Introduction

Raw data

The raw data is shown below. There were 196 rows and 42 columns, consisting of 28 patients sampled over 7 timepoints. However, there is a significant amount of missing data, resulting in only 120 usable datapoints.

missing.ix <- Reduce(intersect, apply(df[,14:42], 2, function(x) which(is.na(x))))
df.raw <- df[-missing.ix,]
df.raw <- df.raw[order(df.raw$PatientID),]
rownames(df.raw) <- 1:nrow(df.raw)
kable(df.raw, 
      digits = 3,
      row.names = T,
      caption = "Raw Data"
) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%", height = "300px") 
Raw Data
PatientID Age AgeGreater60 Sex LowIntermacs InterMACS RVAD Sensitized VAD Indication Device Type Survival Outcome Time num Total PBMC num lymph lymph live lymph CD3 of live lymph CD19 of live lymph CD19+CD27- CD19+CD27+ CD27+38++plasma blasts CD27-38++ transitional CD27-IgD+ mature naive CD27+IgD- switched memory CD27-IgD- switched memory CD27+IgD+ unswitched memory CD27+IgD-IgM+ switched memory CD27+IgD+IgM+ nonswitched memory CD19+27+IgG+IgM- memory CD19+24dim38dim naive mature CD19+24+38++transitional CD19CD24hiCD38-memory CD19+27-38+CD5+transitionals CD19+CD268+ CD268 of +27-38++transitional CD19+CD11b+ CD19+CD5+ CD19+CD27+CD24hi CD19+CD5+CD24hi CD19+CD5+CD11b+ CD19+27+IgD-38++IgG ASC
1 1 65 older Male High 2 No NA BTT HMII alive Alive s/p OHT 0 169154 35496 20.98 99.53 25.37 19.82 86.02 13.98 1.26 2.54 57.81 11.74 28.06 2.39 21.74 14.36 0.72 81.05 2.86 13.60 0.58 95.22 84.83 10.47 4.33 8.50 1.39 3.13 1.09
2 1 65 older Male High 2 No NA BTT HMII alive Alive s/p OHT 1 63082 9915 15.72 99.26 32.45 26.26 90.06 9.94 0.89 2.09 54.27 7.70 35.71 2.32 31.10 20.47 1.57 84.06 2.59 10.79 0.26 97.02 92.59 7.50 3.13 6.42 1.32 2.44 1.49
3 1 65 older Male High 2 No NA BTT HMII alive Alive s/p OHT 3 75921 21721 28.61 99.31 24.86 33.01 86.38 13.62 1.54 1.26 45.74 10.98 40.44 2.84 31.80 18.04 1.15 86.21 4.82 6.56 0.24 90.41 73.33 10.00 7.25 10.12 3.22 5.67 2.05
4 1 65 older Male High 2 No NA BTT HMII alive Alive s/p OHT 8 213808 36002 16.84 99.19 35.95 14.55 67.94 32.06 4.06 2.17 27.92 28.92 39.64 3.52 20.31 9.10 1.33 68.69 7.10 20.42 1.88 87.38 44.25 10.78 11.12 19.15 5.54 7.35 1.07
5 1 65 older Male High 2 No NA BTT HMII alive Alive s/p OHT 21 106970 31867 29.79 99.18 42.36 15.02 84.69 15.31 1.94 2.11 55.41 13.14 29.13 2.32 20.22 12.88 2.08 83.21 5.20 9.50 1.10 94.44 73.00 9.04 6.99 9.22 2.80 4.78 1.44
6 2 81 older Male Low 3 No NA DT HMII dead Died 0 119377 77509 64.93 99.61 58.50 7.99 62.67 37.33 9.10 4.82 36.73 31.90 25.62 5.74 11.47 8.95 0.40 69.54 7.69 9.41 2.61 76.44 57.91 23.29 20.76 15.39 5.42 14.90 2.60
7 2 81 older Male Low 3 No NA DT HMII dead Died 3 96940 57658 59.48 99.70 55.52 9.21 72.15 27.85 3.85 3.82 48.41 21.75 23.49 6.35 12.10 16.04 0.28 75.30 3.72 6.65 2.11 80.52 67.33 15.80 15.23 6.71 1.32 8.86 0.26
8 2 81 older Male Low 3 No NA DT HMII dead Died 5 154373 73969 47.92 99.65 45.18 8.77 67.36 32.64 6.59 2.77 39.10 26.28 27.95 6.67 17.89 15.71 0.39 71.57 6.26 9.88 1.78 81.86 72.07 20.12 20.17 14.80 5.99 14.65 3.71
9 2 81 older Male Low 3 No NA DT HMII dead Died 8 155610 63375 40.73 99.51 38.42 12.40 71.69 28.31 4.36 2.62 37.82 23.43 33.50 5.24 17.37 14.92 0.60 65.04 2.80 11.23 2.54 83.55 37.56 15.55 15.45 9.55 2.66 10.91 0.28
10 2 81 older Male Low 3 No NA DT HMII dead Died 14 244245 115109 47.13 99.47 47.40 9.37 67.50 32.50 3.99 2.08 50.08 26.01 17.22 6.68 15.34 15.25 0.23 77.86 3.21 10.49 1.48 82.13 59.19 21.45 19.77 14.45 2.57 14.94 0.94
11 3 58 younger Male Low 3 Yes No BTT HMII alive Alive s/p OHT 1 246811 86529 35.06 93.49 67.96 17.48 81.66 18.34 2.05 6.85 51.33 13.69 29.91 5.08 15.68 5.94 0.19 61.95 0.76 13.51 4.07 92.62 34.06 6.90 12.07 8.00 2.54 3.85 0.26
12 3 58 younger Male Low 3 Yes No BTT HMII alive Alive s/p OHT 3 406771 90343 22.21 97.34 49.28 17.02 75.04 24.96 4.47 18.36 28.45 21.86 46.07 3.62 24.97 4.33 0.13 47.41 1.64 19.70 18.15 72.76 8.08 25.69 28.35 13.95 4.66 19.20 0.41
13 3 58 younger Male Low 3 Yes No BTT HMII alive Alive s/p OHT 21 330191 81636 24.72 94.57 45.14 14.70 70.51 29.49 8.94 16.23 38.92 23.53 30.83 6.72 23.93 7.07 0.18 46.63 2.35 16.12 8.59 74.45 39.90 24.19 21.18 12.62 3.26 12.93 0.38
14 4 67 older Male Low 3 No Yes BTT HVAD alive Alive s/p OHT 0 309890 41896 13.52 97.03 47.13 8.08 77.31 22.69 2.74 3.29 59.03 19.80 19.37 1.80 21.10 4.46 19.79 71.98 2.41 17.39 0.20 95.80 75.00 5.79 1.89 16.45 0.91 0.67 0.77
15 4 67 older Male Low 3 No Yes BTT HVAD alive Alive s/p OHT 1 274228 15653 5.71 97.85 35.76 28.87 51.45 48.55 2.60 2.13 30.85 43.13 22.77 3.26 20.09 3.92 30.51 57.58 1.31 28.20 0.18 91.23 59.57 12.08 2.04 34.44 0.38 0.16 0.16
16 4 67 older Male Low 3 No Yes BTT HVAD alive Alive s/p OHT 3 591939 57589 9.73 96.76 38.20 42.54 75.71 24.29 2.14 1.33 53.03 21.26 24.16 1.55 25.50 4.02 20.27 67.18 1.14 21.71 0.11 95.28 66.98 3.79 2.43 17.70 0.65 0.84 0.18
17 4 67 older Male Low 3 No Yes BTT HVAD alive Alive s/p OHT 5 189634 46711 24.63 97.74 34.84 37.32 75.30 24.70 0.70 0.55 46.97 22.05 30.28 0.70 23.67 1.56 31.44 73.33 0.75 18.59 0.00 84.61 86.17 3.44 2.11 18.68 0.94 0.63 0.69
18 4 67 older Male Low 3 No Yes BTT HVAD alive Alive s/p OHT 14 182319 40475 22.20 96.28 43.36 11.77 69.34 30.66 1.59 0.89 52.64 27.32 18.56 1.48 22.79 2.88 28.61 75.47 0.33 14.54 0.06 95.90 53.66 7.81 3.27 21.04 1.57 1.50 0.72
19 5 65 older Female High 2 No Yes BTT HMII alive Alive s/p OHT 0 1073919 336538 31.34 100.00 63.29 14.80 93.05 6.95 0.73 8.05 75.49 4.55 19.16 0.80 7.25 8.78 28.04 88.72 2.51 5.18 1.02 90.49 68.18 2.32 4.09 1.69 3.48 1.27 5.78
20 5 65 older Female High 2 No Yes BTT HMII alive Alive s/p OHT 1 270301 97229 35.97 91.63 18.43 61.67 97.68 2.32 0.20 5.85 69.45 1.63 28.19 0.73 27.57 21.43 17.05 84.32 7.26 2.73 0.84 99.44 97.60 0.54 2.64 1.51 3.63 0.29 3.40
21 5 65 older Female High 2 No Yes BTT HMII alive Alive s/p OHT 3 499016 78019 15.63 89.41 25.83 51.38 92.98 7.02 0.32 4.66 51.34 6.31 41.59 0.76 36.49 5.96 22.81 73.48 6.73 7.93 0.33 98.78 96.52 1.31 2.06 5.19 2.68 0.59 2.45
22 5 65 older Female High 2 No Yes BTT HMII alive Alive s/p OHT 8 115825 32513 28.07 91.04 44.61 15.76 80.52 19.48 4.16 7.82 49.91 18.37 30.52 1.20 14.87 3.33 30.19 61.08 8.70 12.56 1.58 89.80 76.71 2.87 6.07 11.12 8.08 1.59 1.17
23 5 65 older Female High 2 No Yes BTT HMII alive Alive s/p OHT 21 407080 97651 23.99 91.29 66.96 12.20 84.66 15.34 2.88 5.60 61.89 13.73 22.65 1.74 20.16 7.23 29.39 68.49 5.28 10.60 1.17 92.50 54.35 3.48 6.53 9.61 8.34 1.85 2.83
24 6 43 younger Male High 2 Yes No BTT HMII alive Alive s/p OHT 0 515760 266115 51.60 96.30 69.66 8.57 72.22 27.78 6.26 4.77 43.09 26.01 28.26 2.63 16.67 2.67 28.25 50.35 0.41 12.78 4.47 84.76 19.47 9.56 4.55 8.29 0.87 2.18 0.75
25 6 43 younger Male High 2 Yes No BTT HMII alive Alive s/p OHT 8 539521 223872 41.49 96.72 60.02 12.63 80.13 19.87 1.66 2.48 43.56 19.52 35.24 1.68 20.46 3.67 20.54 55.76 0.18 10.80 2.43 87.82 13.25 6.55 5.96 4.65 0.95 2.51 0.65
26 6 43 younger Male High 2 Yes No BTT HMII alive Alive s/p OHT 21 564317 216312 38.33 97.40 55.80 4.83 74.19 25.81 5.35 5.77 37.90 25.72 34.76 1.62 17.99 2.18 31.75 61.61 0.77 13.44 6.97 76.55 6.97 10.34 11.14 4.46 0.94 5.14 2.43
27 7 49 younger Male Low 3 No No BTT HVAD alive Alive s/p OHT 0 89251 65181 73.03 76.76 84.53 1.15 74.65 25.35 1.91 1.39 55.03 22.40 19.62 2.95 15.65 6.80 25.85 70.49 0.87 23.26 0.23 83.16 25.00 8.85 5.03 11.63 1.56 2.95 3.08
28 7 49 younger Male Low 3 No No BTT HVAD alive Alive s/p OHT 3 375167 234817 62.59 98.46 60.29 10.18 87.72 12.28 0.46 0.48 68.54 9.90 19.03 2.54 30.58 19.69 15.54 72.77 0.13 20.35 0.03 92.56 36.28 3.75 2.70 9.17 1.19 1.27 0.63
29 7 49 younger Male Low 3 No No BTT HVAD alive Alive s/p OHT 5 256306 162176 63.27 93.82 68.43 8.45 85.45 14.55 0.70 0.53 71.99 11.98 13.37 2.66 32.41 17.96 12.46 69.74 0.41 24.63 0.02 95.72 44.12 4.47 2.23 12.99 1.30 1.20 0.77
30 7 49 younger Male Low 3 No No BTT HVAD alive Alive s/p OHT 8 463136 330714 71.41 98.79 70.51 7.30 87.51 12.49 1.78 0.58 69.60 10.61 17.82 1.97 29.83 14.55 12.59 75.03 0.15 21.61 0.07 93.59 30.94 4.05 3.07 10.00 1.08 1.35 0.62
31 8 43 younger Male Low 3 No No BTT HMII alive Alive s/p OHT 0 92221 27988 30.35 85.45 55.78 1.51 85.56 14.44 0.28 0.28 55.28 4.72 39.72 0.28 0.00 1.89 1.89 89.44 0.00 2.78 0.00 15.00 0.00 3.33 0.83 0.28 0.28 3.33 0.00
32 8 43 younger Male Low 3 No No BTT HMII alive Alive s/p OHT 1 273746 46755 17.08 93.25 50.14 7.57 88.27 11.73 0.09 2.12 54.52 3.36 42.00 0.12 8.40 6.87 0.00 92.00 0.12 3.45 0.00 3.94 11.43 1.12 0.64 2.76 0.33 2.15 0.00
33 8 43 younger Male Low 3 No No BTT HMII alive Alive s/p OHT 3 156053 71752 45.98 78.38 35.43 27.12 78.57 21.43 0.24 0.60 58.56 5.80 18.33 17.31 12.67 72.96 0.53 94.06 0.49 5.11 0.28 22.79 27.47 0.81 5.53 0.01 0.29 0.19 0.00
34 8 43 younger Male Low 3 No No BTT HMII alive Alive s/p OHT 8 156053 71752 45.98 78.38 35.43 27.12 78.57 21.43 0.24 0.60 58.56 5.80 18.33 17.31 12.67 72.96 0.53 94.06 0.49 5.11 0.28 22.79 27.47 0.81 5.53 0.01 0.29 0.19 0.00
35 9 60 younger Female High 2 Yes Yes BTT PVAD dead Died 3 116462 41076 35.27 95.44 33.79 25.24 84.51 15.49 2.92 0.99 45.34 13.19 39.12 2.36 17.56 8.75 19.73 79.35 5.17 10.80 0.71 90.04 29.59 10.16 6.81 8.21 3.37 4.86 4.47
36 9 60 younger Female High 2 Yes Yes BTT PVAD dead Died 5 222650 85649 38.47 98.40 32.09 19.82 77.54 22.46 5.02 1.10 43.96 19.20 33.50 3.34 17.47 9.77 15.63 71.89 6.68 14.08 0.70 86.49 27.87 13.04 9.84 12.30 4.35 7.98 6.39
37 9 60 younger Female High 2 Yes Yes BTT PVAD dead Died 8 102274 63057 61.65 97.73 27.89 20.66 73.75 26.25 7.33 3.45 35.79 21.92 37.73 4.56 9.97 8.20 16.96 63.96 2.40 11.16 3.07 73.97 5.92 16.69 17.77 11.74 4.90 12.19 6.09
38 9 60 younger Female High 2 Yes Yes BTT PVAD dead Died 14 133530 81665 61.16 96.46 22.03 16.15 35.30 64.70 13.74 2.33 8.48 52.28 26.37 12.87 5.29 5.79 19.15 47.48 5.65 17.83 5.43 39.88 6.08 52.08 62.11 37.45 21.95 50.15 3.84
39 10 36 younger Male High 2 Yes NA BTT PVAD alive Alive s/p OHT 0 481928 162052 33.63 99.70 64.73 5.05 3.91 96.09 8.89 1.34 0.72 94.99 3.14 1.15 4.07 0.84 0.32 90.28 2.84 1.58 1.27 6.82 9.17 14.61 90.11 4.95 5.71 13.01 0.51
40 10 36 younger Male High 2 Yes NA BTT PVAD alive Alive s/p OHT 8 38113 21486 56.37 97.86 58.44 7.88 20.47 79.53 8.09 1.69 1.21 75.00 18.96 4.83 10.07 2.59 0.53 82.19 2.29 9.00 3.59 15.16 0.00 23.07 63.95 4.77 6.52 16.91 0.40
41 10 36 younger Male High 2 Yes NA BTT PVAD alive Alive s/p OHT 21 80991 19210 23.72 99.00 50.61 3.50 12.33 87.67 14.74 4.81 1.05 85.41 11.13 2.41 8.48 0.52 1.56 80.00 8.27 3.76 7.41 13.68 0.00 30.83 77.74 10.83 14.14 19.85 2.81
42 11 61 older Male High 2 No No BTT HMII alive Alive s/p OHT 0 411616 119210 28.96 90.82 73.63 5.88 74.30 25.70 2.80 6.46 57.94 19.20 18.72 4.13 19.10 15.42 15.30 53.51 13.78 25.19 0.21 47.23 53.77 5.75 5.92 22.63 3.46 2.55 0.74
43 11 61 older Male High 2 No No BTT HMII alive Alive s/p OHT 1 430509 79492 18.46 89.31 52.10 20.45 87.13 12.87 1.23 4.76 58.11 9.15 30.55 2.20 29.86 18.05 10.37 63.63 10.23 17.43 0.06 74.05 72.50 2.36 2.98 10.68 1.72 1.03 2.89
44 11 61 older Male High 2 No No BTT HMII alive Alive s/p OHT 3 395713 119670 30.24 97.90 23.79 6.27 82.17 17.83 0.79 1.99 49.95 13.01 33.91 3.13 28.36 16.83 5.00 56.33 7.13 26.73 0.02 79.50 82.88 1.96 5.55 16.30 3.61 1.88 0.42
45 11 61 older Male High 2 No No BTT HMII alive Alive s/p OHT 8 329840 73468 22.27 98.63 70.61 3.73 51.50 48.50 8.32 3.77 33.46 41.29 20.92 4.33 25.19 5.89 24.07 46.40 3.73 34.64 0.07 64.14 28.43 10.35 9.61 38.78 4.62 5.55 2.91
46 11 61 older Male High 2 No No BTT HMII alive Alive s/p OHT 14 350294 88823 25.36 97.63 74.37 4.37 73.35 26.65 4.91 4.83 56.54 21.32 19.84 2.30 34.68 11.25 18.37 61.98 5.78 20.55 0.04 63.91 47.54 6.81 6.17 21.77 2.90 2.74 3.00
47 12 38 younger Male High 2 Yes Yes BTT TAH dead Died s/p OHT 0 128451 60927 47.43 96.82 60.77 4.70 77.26 22.74 2.17 1.95 55.58 18.40 21.54 4.47 27.86 18.89 15.33 69.90 0.94 24.40 0.14 70.73 1.85 8.01 8.66 20.21 5.45 3.32 0.39
48 12 38 younger Male High 2 Yes Yes BTT TAH dead Died s/p OHT 1 142796 22573 15.81 80.67 47.31 2.28 63.37 36.63 9.16 6.51 26.02 33.49 37.11 3.37 38.13 7.50 13.13 53.25 5.78 26.27 1.57 40.24 0.00 27.71 22.17 29.16 18.80 15.66 8.11
49 12 38 younger Male High 2 Yes Yes BTT TAH dead Died s/p OHT 3 157167 38008 24.18 93.47 34.86 3.95 69.71 30.29 3.92 3.92 33.64 26.30 35.78 4.28 39.54 11.49 14.94 66.86 2.85 24.31 1.04 64.50 7.27 21.60 20.53 27.66 14.54 14.18 2.12
50 12 38 younger Male High 2 Yes Yes BTT TAH dead Died s/p OHT 5 121823 23097 18.96 93.48 2.51 2.87 99.52 0.48 0.00 0.00 0.97 0.48 98.55 0.00 50.00 25.00 25.00 94.35 0.00 2.58 0.00 0.16 0.00 17.45 87.40 0.00 0.16 18.26 0.00
51 12 38 younger Male High 2 Yes Yes BTT TAH dead Died s/p OHT 8 196515 88579 45.07 85.54 17.17 3.28 59.38 40.62 3.21 1.41 35.44 35.28 23.87 5.42 62.85 12.32 3.71 66.05 4.38 23.22 0.27 74.53 31.43 25.51 28.12 39.78 25.79 21.57 0.79
52 13 61 older Male High 1 Yes No BTT TAH alive Alive s/p OHT 1 87827 15751 17.93 76.38 31.60 12.17 77.46 22.54 6.35 2.53 50.34 20.22 27.05 2.39 19.68 7.62 19.05 75.55 7.38 9.49 0.71 94.26 89.19 11.68 17.96 18.78 6.97 10.86 12.46
53 13 61 older Male High 1 Yes No BTT TAH alive Alive s/p OHT 3 55539 13608 24.50 79.00 28.02 15.78 81.90 18.10 6.13 2.42 66.92 14.33 14.98 3.77 20.81 17.45 18.12 78.01 7.08 8.61 0.72 94.04 70.73 11.32 16.69 15.21 5.78 10.50 16.39
54 13 61 older Male High 1 Yes No BTT TAH alive Alive s/p OHT 8 295005 36785 12.47 94.58 36.55 9.14 86.32 13.68 2.45 1.82 37.01 12.52 49.18 1.29 16.24 4.71 12.47 8.36 1.73 6.64 0.81 84.78 46.55 12.52 9.28 5.94 1.73 4.47 2.49
55 13 61 older Male High 1 Yes No BTT TAH alive Alive s/p OHT 14 60695 20471 33.73 94.91 6.93 7.42 83.97 16.03 3.40 0.69 24.36 14.71 59.33 1.60 18.43 5.07 11.52 54.82 2.91 16.10 0.50 82.44 30.00 51.08 13.25 10.41 3.68 11.31 6.06
56 14 39 younger Female Low 3 Yes No BTT HMII dead Died post OHT 0 52016 30920 59.44 95.77 55.57 24.62 58.73 41.27 2.14 0.97 53.01 5.05 5.45 36.50 9.39 81.29 2.88 81.48 5.65 9.02 0.48 97.75 71.83 0.67 22.89 0.34 1.78 0.49 NA
57 14 39 younger Female Low 3 Yes No BTT HMII dead Died post OHT 1 75516 48801 64.62 92.13 52.90 18.64 70.39 29.61 2.92 2.21 62.18 5.13 7.98 24.70 13.79 69.54 3.19 84.14 3.64 6.36 0.41 94.93 43.78 0.57 15.61 0.23 0.91 0.36 NA
58 14 39 younger Female Low 3 Yes No BTT HMII dead Died post OHT 3 250612 130454 52.05 95.75 60.05 3.29 46.04 53.96 1.34 3.89 33.08 32.47 8.29 26.15 3.56 34.61 0.09 76.71 0.39 13.20 0.28 39.79 14.38 2.26 1.56 10.45 2.89 3.04 0.08
59 14 39 younger Female Low 3 Yes No BTT HMII dead Died post OHT 8 212478 120477 56.70 85.12 48.36 12.51 79.57 20.43 1.35 2.14 58.73 12.87 18.44 9.96 5.63 34.32 0.11 86.74 0.41 6.88 0.15 35.30 53.09 1.23 0.75 2.36 1.04 1.47 NA
60 15 70 older Female High 2 No NA DT HMII alive Alive 0 86135 70306 81.62 93.71 74.50 1.59 61.74 38.26 1.91 2.67 33.97 14.50 26.81 24.71 20.67 62.47 9.03 73.00 2.00 24.71 1.27 20.23 71.43 0.48 12.79 0.38 3.24 0.48 4.60
61 15 70 older Female High 2 No NA DT HMII alive Alive 1 81004 41334 51.03 94.06 53.64 2.99 61.70 38.30 3.10 2.75 39.07 13.34 22.12 25.47 12.58 63.11 15.99 78.92 2.93 17.38 3.35 23.24 71.88 0.52 11.70 0.09 1.29 0.52 9.83
62 15 70 older Female High 2 No NA DT HMII alive Alive 5 90024 39434 43.80 95.32 68.83 5.62 71.59 28.41 0.85 1.66 35.23 13.02 35.80 15.96 22.94 52.37 5.85 76.42 1.14 21.83 0.62 14.54 48.57 0.38 4.69 0.00 0.80 0.38 0.62
63 15 70 older Female High 2 No NA DT HMII alive Alive 14 106791 80532 75.41 88.82 79.92 1.57 62.32 37.68 1.70 4.02 28.04 22.86 34.11 15.00 23.92 39.41 15.03 61.61 3.66 34.55 1.97 20.36 73.33 0.80 8.57 0.27 1.79 0.98 0.71
64 16 63 older Male High 1 No NA BTT HMII dead Died 3 366713 77282 21.07 98.61 66.55 12.44 93.65 6.35 0.53 0.60 62.22 4.68 32.38 0.72 19.94 9.89 28.23 83.18 0.51 14.38 0.30 94.72 54.39 2.07 11.33 2.86 3.40 1.48 3.58
65 16 63 older Male High 1 No NA BTT HMII dead Died 5 262251 55344 21.10 99.55 62.37 14.39 93.00 7.00 1.03 0.69 65.71 4.74 28.37 1.17 25.93 17.28 25.93 79.87 1.20 17.05 0.48 97.76 58.18 2.27 11.13 4.12 4.91 1.87 6.10
66 16 63 older Male High 1 No NA BTT HMII dead Died 8 309193 61091 19.76 98.15 56.62 7.24 83.61 16.39 3.27 1.52 58.03 11.90 28.13 1.93 24.65 11.36 1.52 79.53 1.45 13.33 0.97 92.63 25.76 3.50 9.23 7.57 3.73 3.48 0.00
67 17 68 older Male Low 3 No NA DT HMII dead Died 0 147997 54649 36.93 91.75 79.63 7.62 74.86 25.14 0.60 3.45 62.45 22.87 12.30 2.38 45.98 9.28 13.61 74.33 2.67 21.19 0.84 91.31 90.15 3.58 2.25 21.79 1.36 0.58 0.34
68 17 68 older Male Low 3 No NA DT HMII dead Died 1 177687 56400 31.74 97.13 69.76 12.60 81.83 18.17 0.65 1.61 55.77 16.13 26.02 2.09 53.19 10.70 11.57 82.92 1.38 13.98 0.28 89.83 85.59 3.78 1.85 14.65 1.30 0.68 0.18
69 17 68 older Male Low 3 No NA DT HMII dead Died 3 233155 67140 28.80 98.44 74.89 5.32 81.96 18.04 1.00 1.88 49.03 16.59 32.90 1.48 49.92 7.20 12.83 83.44 1.39 12.98 0.14 89.70 69.70 4.58 1.62 13.89 1.08 1.08 0.17
70 17 68 older Male Low 3 No NA DT HMII dead Died 5 99257 29455 29.68 97.33 62.48 4.73 80.60 19.40 1.84 1.25 45.06 18.36 35.25 1.33 44.61 6.32 12.27 82.89 0.88 13.20 0.09 83.55 35.29 7.23 3.54 14.75 2.29 3.02 1.20
71 17 68 older Male Low 3 No NA DT HMII dead Died 8 76358 18702 24.49 97.32 29.48 11.10 86.29 13.71 2.43 0.50 37.92 12.13 48.32 1.63 46.45 11.70 11.35 83.96 0.20 10.50 0.06 57.82 50.00 9.95 29.60 10.20 0.64 9.36 0.80
72 18 37 younger Female High 2 Yes NA BTT PVAD dead Died s/p OHT 0 56374 12141 21.54 67.38 27.12 3.32 51.84 48.16 1.84 1.84 2.94 45.59 47.79 3.68 0.00 0.00 2.29 92.65 0.37 5.51 0.72 30.15 60.00 6.99 9.56 0.37 0.00 3.68 0.00
73 18 37 younger Female High 2 Yes NA BTT PVAD dead Died s/p OHT 3 7018 3464 49.36 89.03 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
74 18 37 younger Female High 2 Yes NA BTT PVAD dead Died s/p OHT 8 80210 21828 27.21 85.74 20.38 3.50 61.98 38.02 8.55 2.90 5.50 34.96 56.18 3.36 0.00 0.00 3.21 89.31 0.61 2.60 1.49 23.66 0.00 21.37 9.16 1.07 0.15 5.34 1.31
75 18 37 younger Female High 2 Yes NA BTT PVAD dead Died s/p OHT 14 55068 19145 34.77 75.16 16.90 3.30 64.00 36.00 4.21 2.11 10.53 34.11 52.42 2.95 0.00 0.00 2.94 92.42 0.00 3.37 1.00 22.74 0.00 12.42 11.16 0.63 0.00 4.63 0.00
76 18 37 younger Female High 2 Yes NA BTT PVAD dead Died s/p OHT 21 76673 23615 30.80 83.39 26.27 2.06 63.55 36.45 5.42 5.42 5.17 34.48 58.13 2.22 0.00 0.00 2.72 89.41 0.00 3.45 1.95 18.23 4.55 18.23 8.87 1.23 0.25 4.68 0.00
77 19 46 younger Male High 2 No NA DT HMII alive Alive 0 226762 25868 11.41 89.41 25.31 12.05 55.06 44.94 11.31 2.69 14.79 42.57 39.95 2.69 19.87 3.88 32.88 68.95 10.48 13.57 0.59 80.62 22.67 31.91 19.38 15.58 10.37 16.55 8.18
78 19 46 younger Male High 2 No NA DT HMII alive Alive 3 244049 39386 16.14 94.94 40.46 13.23 40.93 59.07 13.36 3.17 8.33 54.90 32.28 4.49 18.52 4.21 23.83 55.18 13.87 19.39 1.44 70.53 35.67 33.62 33.47 22.86 15.26 25.81 5.11
79 19 46 younger Male High 2 No NA DT HMII alive Alive 5 344736 250121 72.55 89.83 7.54 3.81 45.84 54.16 3.96 1.53 13.39 52.70 32.18 1.72 67.14 2.92 2.44 74.68 6.81 13.91 0.21 63.33 58.02 33.05 30.94 19.83 13.96 26.74 0.22
80 19 46 younger Male High 2 No NA DT HMII alive Alive 8 216934 140093 64.58 96.37 4.78 7.22 22.49 77.51 13.49 0.64 7.22 67.61 15.17 10.00 35.88 7.82 17.22 51.28 39.40 6.56 1.33 62.59 58.06 66.67 64.59 53.00 50.26 63.78 4.97
81 19 46 younger Male High 2 No NA DT HMII alive Alive 14 215820 66785 30.94 92.54 17.20 6.01 25.88 74.12 9.15 1.45 3.09 71.72 22.44 2.74 65.43 3.42 0.04 69.38 15.93 8.85 0.74 50.58 24.07 52.78 48.88 32.28 28.19 44.04 0.04
82 19 46 younger Male High 2 No NA DT HMII alive Alive 21 670692 523014 77.98 96.05 5.07 2.64 23.61 76.39 4.11 0.49 4.79 73.53 18.68 3.00 42.40 2.72 4.49 71.95 20.81 5.51 0.45 48.22 44.62 60.13 60.26 39.07 36.82 56.38 1.49
83 20 75 older Male High 2 No NA DT HMII alive Alive 3 121794 28129 23.10 97.05 52.79 16.42 70.08 29.92 0.29 0.54 50.71 21.46 22.96 4.86 34.18 12.51 23.49 39.78 0.31 55.18 0.35 96.01 16.67 1.74 5.27 23.38 2.95 1.00 0.10
84 20 75 older Male High 2 No NA DT HMII alive Alive 5 82157 17295 21.05 99.20 72.54 7.41 68.21 31.79 1.10 0.55 48.47 22.90 23.60 5.04 31.90 13.57 18.10 37.84 0.47 56.73 0.35 92.37 14.29 2.05 7.79 26.20 4.72 1.42 1.02
85 20 75 older Male High 2 No NA DT HMII alive Alive 8 154810 27067 17.48 98.04 72.18 5.03 74.91 25.09 0.90 1.65 58.05 14.76 20.82 6.37 23.43 25.71 22.57 44.94 0.60 49.36 0.71 86.59 22.73 3.30 7.64 20.15 4.87 1.72 0.00
86 20 75 older Male High 2 No NA DT HMII alive Alive 21 197936 15446 7.80 95.86 50.35 5.44 79.13 20.87 3.98 3.11 52.30 15.90 29.81 1.99 20.23 5.78 24.86 40.62 1.12 50.06 1.89 62.36 8.00 7.20 9.94 13.17 3.60 5.09 1.55
87 21 30 younger Female Low 3 No Yes BTT HMII alive Alive s/p OHT 0 114586 35996 31.41 100.00 59.03 11.31 75.87 24.13 1.52 7.69 69.83 11.67 10.20 8.30 27.80 29.94 2.04 84.64 0.64 10.54 0.61 91.57 61.34 4.99 7.15 11.57 5.92 2.53 1.92
88 21 30 younger Female Low 3 No Yes BTT HMII alive Alive s/p OHT 3 349054 129419 37.08 95.25 62.35 8.81 79.64 20.36 2.72 4.45 57.70 13.93 25.89 2.48 57.14 13.49 8.69 74.29 3.76 12.35 1.41 85.87 62.53 3.41 5.43 7.67 1.40 1.92 1.06
89 21 30 younger Female Low 3 No Yes BTT HMII alive Alive s/p OHT 5 447309 176789 39.52 93.00 53.57 9.32 80.92 19.08 2.16 5.73 56.72 12.51 28.43 2.34 47.68 12.19 11.00 73.80 4.44 12.53 2.57 83.88 49.09 5.24 6.82 8.32 2.18 2.28 1.72
90 21 30 younger Female Low 3 No Yes BTT HMII alive Alive s/p OHT 8 558290 301468 54.00 95.62 47.73 7.04 79.13 20.87 3.43 7.49 53.19 13.48 30.32 3.01 36.36 14.54 12.10 70.55 5.15 11.89 3.51 79.79 52.89 6.98 8.50 7.82 1.92 3.09 1.96
91 21 30 younger Female Low 3 No Yes BTT HMII alive Alive s/p OHT 14 479402 362845 75.69 97.95 44.30 6.30 76.53 23.47 3.90 11.10 34.99 12.04 50.31 2.66 21.93 14.88 10.05 70.21 4.66 7.31 7.73 49.96 29.38 13.44 17.37 5.64 3.01 6.95 2.28
92 21 30 younger Female Low 3 No Yes BTT HMII alive Alive s/p OHT 21 334942 191418 57.15 93.28 70.35 4.04 70.89 29.11 6.20 17.70 56.45 19.57 20.44 3.53 38.91 14.64 10.35 55.11 6.61 14.34 8.70 80.96 58.92 7.58 7.96 10.11 2.08 2.08 2.48
93 22 67 older Male High 2 Yes NA BTT HMII dead Died 0 129980 41600 32.00 91.40 39.03 2.85 49.63 50.37 7.02 11.55 28.37 29.11 20.79 21.72 19.43 28.34 9.80 69.59 4.90 20.33 6.13 31.15 33.60 20.15 22.27 0.28 4.25 16.54 2.09
94 22 67 older Male High 2 Yes NA BTT HMII dead Died 1 59478 20182 33.93 89.07 44.50 10.08 60.93 39.07 2.10 2.70 42.88 13.69 17.66 25.77 13.06 46.12 3.13 81.35 1.16 16.11 0.83 37.58 6.12 5.52 11.64 0.06 1.49 2.76 4.07
95 22 67 older Male High 2 Yes NA BTT HMII dead Died 5 26646 9754 36.61 98.69 58.32 7.57 64.88 35.12 2.06 0.55 46.50 12.62 17.70 23.18 15.07 53.31 12.50 83.40 1.10 15.64 1.76 28.12 0.00 7.54 16.46 0.55 1.92 6.04 13.00
96 22 67 older Male High 2 Yes NA BTT HMII dead Died 8 73641 21944 29.80 96.79 48.91 4.47 62.74 37.26 2.84 3.16 37.47 17.79 24.95 19.79 19.46 38.11 9.46 79.79 1.89 13.58 1.24 27.05 20.00 14.11 13.37 0.21 2.74 10.53 3.31
97 23 68 older Male High 2 No No BTT HMII alive Alive s/p OHT 0 44004 20612 46.84 43.65 15.61 21.54 49.64 50.36 10.37 0.57 15.94 45.15 33.18 5.73 24.69 5.68 20.25 44.58 1.55 36.17 0.88 67.23 27.27 32.35 15.48 20.02 4.49 10.94 2.76
98 23 68 older Male High 2 No No BTT HMII alive Alive s/p OHT 1 99731 23424 23.49 91.81 41.33 10.89 35.24 64.76 11.83 1.62 18.71 59.21 16.19 5.89 28.16 5.89 18.53 39.98 4.66 42.16 2.33 84.07 2.63 22.55 15.21 40.28 6.07 9.70 3.96
99 23 68 older Male High 2 No No BTT HMII alive Alive s/p OHT 8 129637 50001 38.57 90.65 44.41 6.05 30.66 69.34 17.12 4.05 10.62 66.50 19.60 3.28 26.41 3.46 13.37 39.74 7.15 39.49 11.91 71.97 2.70 32.48 33.14 43.91 7.81 22.70 3.82
100 23 68 older Male High 2 No No BTT HMII alive Alive s/p OHT 14 134368 46363 34.50 88.32 34.25 4.63 32.47 67.53 12.86 8.38 8.91 64.89 22.88 3.32 24.78 3.89 13.58 36.16 6.01 40.33 13.77 72.11 1.89 29.89 26.36 42.80 6.17 15.97 3.29
101 24 64 older Male High 2 No Yes BTT HMII alive Alive s/p OHT 0 59881 8971 14.98 97.40 50.06 9.57 75.84 24.16 3.59 8.13 64.83 13.40 15.67 6.10 19.31 18.81 16.83 76.67 6.10 6.82 2.00 88.16 47.06 10.53 11.24 10.41 9.93 6.94 11.93
102 24 64 older Male High 2 No Yes BTT HMII alive Alive s/p OHT 1 60193 23874 39.66 96.72 58.83 16.24 80.37 19.63 1.87 3.31 68.58 10.62 14.72 6.08 21.33 22.83 3.80 80.61 3.63 8.03 0.33 94.16 57.26 2.96 5.01 10.19 4.72 1.41 1.24
103 24 64 older Male High 2 No Yes BTT HMII alive Alive s/p OHT 5 135335 27834 20.57 97.35 44.75 5.61 63.22 36.78 8.49 17.63 39.01 23.16 30.72 7.11 17.53 10.38 12.16 42.70 4.93 19.28 4.68 67.37 4.10 26.64 17.76 19.21 14.54 14.47 10.60
104 24 64 older Male High 2 No Yes BTT HMII alive Alive s/p OHT 8 247113 62042 25.11 96.41 43.93 14.82 78.42 21.58 2.52 5.24 65.35 12.22 17.53 4.90 20.39 17.20 2.93 77.86 2.11 8.78 2.94 84.64 9.70 10.64 13.93 9.90 11.59 8.48 3.51
105 24 64 older Male High 2 No Yes BTT HMII alive Alive s/p OHT 14 142053 56983 40.11 93.01 62.72 12.87 85.10 14.90 0.69 1.64 68.02 8.40 21.03 2.55 17.03 13.68 1.08 83.33 0.37 7.38 0.26 79.00 16.07 3.01 5.85 5.16 5.24 1.77 0.70
106 24 64 older Male High 2 No Yes BTT HMII alive Alive s/p OHT 21 60384 17029 28.20 93.65 40.96 21.36 83.80 16.20 1.47 2.32 63.52 9.69 23.80 2.99 23.55 15.58 0.91 78.16 1.09 8.13 1.00 83.65 17.72 4.70 8.42 5.99 7.22 3.76 1.56
107 25 25 younger Female Low 3 No Yes BTT HMII alive Alive s/p OHT 0 47645 37377 78.45 94.61 2.47 9.41 78.56 17.02 1.71 4.63 34.94 16.09 44.77 4.21 4.35 28.99 59.42 67.26 3.10 26.22 7.50 72.82 53.90 13.50 28.20 3.04 7.70 7.67 0.96
108 25 25 younger Female Low 3 No Yes BTT HMII alive Alive s/p OHT 1 87452 43025 49.20 91.54 1.98 26.71 92.06 6.76 0.77 4.67 47.21 4.70 45.17 2.92 6.06 18.18 60.61 81.01 3.95 14.10 12.04 85.97 79.84 5.82 14.15 1.11 4.77 2.84 0.00
109 25 25 younger Female Low 3 No Yes BTT HMII alive Alive s/p OHT 3 124056 80873 65.19 95.41 0.37 14.36 90.47 5.98 0.39 0.60 41.84 4.26 51.08 2.81 18.82 35.29 43.53 86.97 0.10 11.32 1.10 75.49 12.12 3.80 14.94 0.24 0.88 3.25 0.00
110 25 25 younger Female Low 3 No Yes BTT HMII alive Alive s/p OHT 8 87840 57990 66.02 95.60 0.11 15.80 89.00 9.49 1.56 5.02 45.89 7.18 43.54 3.39 6.10 19.51 56.71 84.23 3.14 9.95 10.70 72.13 52.73 7.06 17.53 0.74 3.13 5.13 0.13
111 26 49 younger Male Low 3 No NA DT HMII alive Alive 0 196315 122053 62.17 95.37 60.66 22.66 91.10 8.90 0.45 4.95 85.57 6.30 5.42 2.71 13.13 28.32 20.05 86.32 2.45 5.42 0.96 98.68 94.19 2.50 2.96 4.63 1.24 0.33 0.89
112 26 49 younger Male Low 3 No NA DT HMII alive Alive 1 91909 37622 40.93 96.10 43.22 38.85 95.48 4.52 0.14 2.19 85.07 3.09 10.33 1.51 21.03 27.23 9.08 88.84 1.13 3.91 0.23 98.77 88.27 1.16 1.16 2.34 0.38 0.21 0.45
113 26 49 younger Male Low 3 No NA DT HMII alive Alive 5 153762 86281 56.11 95.83 65.54 24.24 91.12 8.88 0.65 2.84 81.43 6.30 9.62 2.65 21.14 28.73 15.89 83.00 1.97 6.74 0.36 98.22 93.16 3.20 2.29 5.47 1.11 0.77 1.56
114 26 49 younger Male Low 3 No NA DT HMII alive Alive 8 82406 47485 57.62 92.71 65.75 17.46 91.79 8.21 1.08 1.99 76.35 5.70 15.41 2.54 20.56 26.95 11.99 78.67 1.35 7.97 0.28 97.55 85.62 2.65 2.56 5.19 1.18 0.86 1.79
115 27 62 older Male High 2 No NA BTT CMAG dead Died 0 105070 38521 36.66 98.39 47.84 12.88 60.82 39.18 11.08 2.89 33.87 29.43 26.68 10.01 19.72 21.92 2.99 59.35 14.64 20.23 3.45 81.40 46.10 24.47 24.27 24.04 11.00 17.02 1.51
116 27 62 older Male High 2 No NA BTT CMAG dead Died 1 162336 73080 45.02 99.08 49.35 13.95 74.24 25.76 7.77 3.84 45.21 19.33 28.85 6.61 19.16 23.45 0.81 65.34 14.40 13.50 3.17 86.68 71.39 23.26 17.70 15.21 6.93 13.29 0.88
117 27 62 older Male High 2 No NA BTT CMAG dead Died 5 111264 60487 54.36 99.15 45.94 16.13 73.16 26.84 4.64 1.82 52.04 20.07 20.94 6.95 22.93 23.47 0.50 77.11 9.30 9.58 1.20 88.92 81.82 22.20 19.01 16.98 6.89 13.57 1.36
118 27 62 older Male High 2 No NA BTT CMAG dead Died 8 165870 67935 40.96 98.88 43.01 15.86 70.23 29.77 7.85 1.88 40.39 20.49 29.66 9.46 22.70 29.00 1.88 68.61 9.60 16.11 1.47 89.57 69.00 20.12 16.59 17.10 7.02 12.40 2.06
119 28 72 older Male Low 4 Yes NA BTT HMII dead Died 3 424375 21971 5.18 71.78 56.98 18.53 79.84 20.16 3.08 0.51 57.91 14.07 24.06 3.97 31.74 15.95 13.65 57.94 0.65 28.51 0.17 93.57 53.33 2.87 9.69 15.61 3.87 2.16 3.99
120 28 72 older Male Low 4 Yes NA BTT HMII dead Died 8 177249 103442 58.36 76.40 15.43 5.64 65.13 34.87 5.21 0.20 35.56 27.91 31.90 4.62 26.60 8.28 13.24 67.89 0.29 20.13 0.00 78.60 22.22 5.50 16.09 25.56 6.80 4.73 0.89

Identification of clusters

We double standardized the data and biclustered it. We found 3 clusters of patients, and 3 clusters of B-cells.

require(pheatmap, quietly = T)
require(RColorBrewer, quietly = T)


colorfun <- function(nlevels, ...){
    if(nlevels > 10) return(standardColors(nlevels))
    if(nlevels <=10) return(pal_d3(...)(nlevels))
} 

annotation.col <- df.raw[,c(1,3,4,6,7,8,9,10,11,12)]
annotation.colors <- lapply(colnames(annotation.col), function(nam){
    colrs <- colorfun(nlevels(annotation.col[[nam]]))
    names(colrs) <- levels(annotation.col[[nam]])
    return(colrs)
})
names(annotation.colors) <- colnames(annotation.col)

df.patient <- double_standardize(t(df.raw[,c(2,14:42)]))
pheatmap(df.patient, 
         color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdBu")))(100),
         breaks = c(seq(-max(abs(df.patient), na.rm=T), 0, length.out = 50), seq(0.001, max(abs(df.patient), na.rm=T), length.out = 50)),
         #scale = "row",
         fontsize = 7,
         cutree_rows = 3,
         cutree_cols = 3,
         annotation_col = annotation.col,
         annotation_colors = annotation.colors
)

Metadata associations

We computed the statistical associations between sample groups specified by the metadata.

Survival

We computed metadata factors that were statistically associated with survival.

df.factors <- df[match(levels(df$PatientID), df$PatientID),]
isfactor <- which(sapply(df.factors, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Survival","Outcome", "LowIntermacs"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$Survival)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)
OR lowerCL upperCL pvalue qvalue Factor
1 4.900 0.7440 39.00 0.0974 0.681 RVAD
2 NA NA NA 0.354 0.823 Device Type
3 2.610 0.1070 188.00 0.569 0.823 Sensitized
4 NA NA NA 0.616 0.823 InterMACS
5 0.677 0.0853 5.94 0.674 0.823 Sex
6 1.480 0.2450 9.82 0.705 0.823 AgeGreater60
7 0.879 0.0656 7.89 1 1 VAD Indication

InterMACS

We computed metadata factors that were statistically associated with a binarized InterMACS score.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("LowIntermacs", "InterMACS"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$LowIntermacs)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)
OR lowerCL upperCL pvalue qvalue Factor
1 NA NA NA 0.135 0.983 Device Type
2 3.070 0.520 21.10 0.246 0.983 AgeGreater60
3 0.583 0.062 5.43 0.653 1 VAD Indication
4 NA NA NA 0.689 1 Outcome
5 1.830 0.288 14.60 0.689 1 RVAD
6 1.310 0.116 15.90 1 1 Sensitized
7 1.210 0.140 9.40 1 1 Sex
8 0.956 0.153 6.40 1 1 Survival

Sensitization

We computed metadata factors that were statistically associated with sensitization.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Sensitized"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$Sensitized)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)
OR lowerCL upperCL pvalue qvalue Factor
1 0.127 0.00192 1.99 0.119 1 Sex
2 2.610 0.10700 188.00 0.569 1 Survival
3 NA NA NA 0.569 1 Outcome
4 0.426 0.02530 5.03 0.608 1 RVAD
5 NA NA NA 0.843 1 Device Type
6 1.310 0.11600 15.90 1 1 LowIntermacs
7 1.230 0.10100 15.40 1 1 AgeGreater60
8 NA NA NA 1 1 InterMACS
9 0.000 0.00000 Inf 1 1 VAD Indication

Sex

We computed metadata factors that were statistically associated with sex.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Sex"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$Sex)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)
OR lowerCL upperCL pvalue qvalue Factor
1 0.127 0.00192 1.99 0.119 0.89 Sensitized
2 3.860 0.48600 49.80 0.198 0.89 AgeGreater60
3 NA NA NA 0.368 1 Outcome
4 NA NA NA 0.447 1 Device Type
5 0.677 0.08530 5.94 0.674 1 RVAD
6 0.677 0.08530 5.94 0.674 1 Survival
7 1.210 0.14000 9.40 1 1 LowIntermacs
8 NA NA NA 1 1 InterMACS
9 1.840 0.15200 103.00 1 1 VAD Indication

Age

We computed metadata factors that were statistically associated with age.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("AgeGreater60"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$AgeGreater60)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)
OR lowerCL upperCL pvalue qvalue Factor
1 0.228 0.0274 1.45 0.114 0.442 RVAD
2 NA NA NA 0.128 0.442 InterMACS
3 NA NA NA 0.159 0.442 Outcome
4 3.860 0.4860 49.80 0.198 0.442 Sex
5 3.070 0.5200 21.10 0.246 0.442 LowIntermacs
6 NA NA NA 0.35 0.525 Device Type
7 1.950 0.2230 25.90 0.655 0.794 VAD Indication
8 1.480 0.2450 9.82 0.705 0.794 Survival
9 1.230 0.1010 15.40 1 1 Sensitized

VAD Indication

We computed metadata factors that were statistically associated with VAD Indication.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("VAD Indication"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$`VAD Indication`)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = F, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)
OR lowerCL upperCL pvalue qvalue Factor
1 NA NA NA 0.000242 0.00217 Outcome
2 0.000 0.0000 1.34 0.0619 0.278 RVAD
3 0.583 0.0620 5.43 0.653 1 LowIntermacs
4 1.950 0.2230 25.90 0.655 1 AgeGreater60
5 NA NA NA 0.837 1 InterMACS
6 NA NA NA 0.877 1 Device Type
7 1.840 0.1520 103.00 1 1 Sex
8 0.000 0.0000 Inf 1 1 Sensitized
9 0.879 0.0656 7.89 1 1 Survival

RVAD

We computed metadata factors that were statistically associated with RVAD.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("RVAD"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$RVAD)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)
OR lowerCL upperCL pvalue qvalue Factor
1 NA NA NA 0.00868 0.0782 Device Type
2 0.000 0.0000 1.34 0.0619 0.205 VAD Indication
3 NA NA NA 0.0739 0.205 Outcome
4 4.900 0.7440 39.00 0.0974 0.205 Survival
5 0.228 0.0274 1.45 0.114 0.205 AgeGreater60
6 NA NA NA 0.327 0.49 InterMACS
7 0.426 0.0253 5.03 0.608 0.689 Sensitized
8 0.677 0.0853 5.94 0.674 0.689 Sex
9 1.830 0.2880 14.60 0.689 0.689 LowIntermacs

Device-type

We computed metadata factors that were statistically associated with Device Type.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[!(names(isfactor) %in% c("Device Type"))], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$`Device Type`)
                      })

fisher.p <- sapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = T, B = 10000)
    ft$p
})

fisher.OR.0 <- lapply(contingency, function(this_table){
    ft <- fisher.test(this_table, simulate.p.value = T, B = 10000)
    if(is.null(ft$estimate)) {
        ft$estimate <- NA
        ft$conf.int <- c(NA,NA)
    }
    ret <- list(OR = ft$estimate, lowerCL = ft$conf.int[1], upperCL = ft$conf.int[2])
    return(ret)
})
fisher.OR <- as.data.frame(do.call(rbind, fisher.OR.0))
fisher.OR$pvalue <- fisher.p
fisher.OR$qvalue <- p.adjust(fisher.p, method = "BH")
fisher.all <- as.data.frame(sapply(fisher.OR, as.numeric))
rownames(fisher.all) <- rownames(fisher.OR)

qtable <- signif(fisher.all[order(fisher.all$pvalue),],3)

qtable %>%
    mutate(
        `Factor` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>% 
    kable(escape = F, row.names = T) %>%
    kable_styling(bootstrap_options = c("striped", 
                                        "hover", 
                                        "condensed",
                                        "responsive"),
                  font_size = 12)
OR lowerCL upperCL pvalue qvalue Factor
1 NA NA NA 0.0089 0.0801 RVAD
2 NA NA NA 0.139 0.528 LowIntermacs
3 NA NA NA 0.274 0.528 InterMACS
4 NA NA NA 0.347 0.528 Outcome
5 NA NA NA 0.352 0.528 AgeGreater60
6 NA NA NA 0.352 0.528 Survival
7 NA NA NA 0.443 0.569 Sex
8 NA NA NA 0.835 0.874 Sensitized
9 NA NA NA 0.874 0.874 VAD Indication

Visualization of metadata associations

We used the GGally package to visualize pairwise dependencies between sample groups, separating out survivors and non-survivors (indicated by color).

df.temp <- df.factors
colnames(df.temp) <- make.names(colnames(df), unique = T)
#invisible(
suppressWarnings(
    suppressMessages(
        ggpairs(df.temp, 
                #mapping = aes(color = Survival), 
                columns = colnames(df.temp)[c(2,3,4,5,6,7,8,9,10,11,12)]) 
    ))

#)

Variability of MCS devices

PCA

Using Principal Component Analysis (PCA), we saw large variability between device types, compared to the variability within device types. We also saw large variability between individual patients. None of the other features were clearly separable.

suppressMessages(require(ggbiplot, quietly = T))
require(ggsci, quietly = T)

isna <- unique(unlist(apply(df[,14:42], 2, function(x) which(is.na(x)))))
pca <- prcomp(double_standardize(df[-isna, 14:42]), center = TRUE, scale. = TRUE)

colorfun <- function(grouping, ...){
    if(nlevels(factor(grouping)) > 10) scale_color_discrete(...)
    if(nlevels(factor(grouping)) <=10) scale_color_d3(...)
} 

plotvars <- names(df)[c(10,1:9,11:17)]
plots.pca <- mclapply(plotvars, function(this_var){
    this_groups <- df[-isna, this_var]
    if(is.numeric(this_groups)) this_groups <- cut(this_groups, 4)
    
    ggbiplot(pca,
             groups = this_groups,
             ellipse = TRUE,
             alpha = 0.3,
             varname.size = 1.2) + 
        colorfun(this_groups, name = this_var) +
        ggtitle(paste0(this_var, " variability")) +
        theme_classic() 
}, mc.cores = detectCores()-1)
names(plots.pca) <- plotvars
#plots.pca$PatientID
#plots.pca$`Device Type`

for(ii in 1:length(plots.pca)){
    cat("  \n###", names(plots.pca)[ii], "\n")
    suppressWarnings(print(plots.pca[[ii]]))
    cat("  \n")
}

Device Type

PatientID

Age

AgeGreater60

Sex

LowIntermacs

InterMACS

RVAD

Sensitized

VAD Indication

Survival

Outcome

Time

num Total PBMC

num lymph

lymph

live lymph

Fisher’s exact test

We analyzed the dependency of device types on the other discrete variables using Fisher’s exact test.

isfactor <- which(sapply(df, is.factor))[-1]
contingency <- lapply(isfactor[names(isfactor) != "Device Type"], 
                      function(this_factor){
                          table(df.factors[,this_factor], df.factors$`Device Type`)
                      })

fisher.p <- sapply(contingency, function(this_table){
    fisher.test(this_table, simulate.p.value = T, B = 10000)$p
})
fisher.q <- p.adjust(fisher.p, method = "BH")

signif.ix <- which(fisher.p < 0.05)
signif.order <- sort(fisher.q[signif.ix], index.return = T)$ix
for(this_ix in signif.ix[signif.order]){
    cat("  \n###", names(contingency)[this_ix], "\n")
    cat(paste0("Benjamini-Hochberg qvalue = ", signif(fisher.q[this_ix], 2)),
        ". \n")
    print(kable(contingency[this_ix][[1]],  row.names = T) %>%
              kable_styling(bootstrap_options = c("striped", 
                                                  "hover", 
                                                  "condensed",
                                                  "responsive"),
                            font_size = 12)
    )
    cat("  \n")
}

RVAD

Benjamini-Hochberg qvalue = 0.077 .
HMII CMAG HVAD PVAD TAH
No 15 1 2 0 0
Yes 5 0 0 3 2

One-way repeated measures anova

We first analyzed the differences in B-cell levels across device types using a one-way repeated measures anova. Here we report any variables that had a statistically significant variance (\(p<0.05\)) across devices, time, or their interaction.

suppressMessages(require(lmerTest, quietly = T))
suppressMessages(require(car, quietly = TRUE))

df.device <- df
colnames(df.device) <- make.names(colnames(df), unique = T)

varnames <- colnames(df.device)[14:42]
models.device <- mclapply(varnames, function(this_var){
    this_formula <- as.formula(paste0(this_var, " ~ Device.Type + (1|PatientID)"))
    invisible(suppressMessages(this_model <- lmer(this_formula, data = df.device)))
    this_anova <- Anova(this_model, type = 2)
    pvals <- this_anova$`Pr(>Chisq)`[c(1)]
    return(list(model = this_model,
                pvals = pvals))
}, mc.cores = detectCores()-1)
names(models.device) <- colnames(df)[14:42]

pvals <- do.call(rbind, lapply(models.device, function(x) x$pvals))
rownames(pvals) <- colnames(df)[14:42]
colnames(pvals) <- c("pvalue")

qBH <- matrix(p.adjust(pvals, method = "BH"), nrow = nrow(pvals))
rownames(qBH) <- colnames(df)[14:42]
colnames(qBH) <- c("qvalue")

sigvars <- apply(apply(pvals, 2, function(x) x<=0.05), 1, function(x) sum(x) > 0)

sig.qtable <- cbind(pvals,qBH)[sigvars,,drop=F][order(apply(pvals[sigvars,,drop=F], 1, min)),,drop=F]

qtable <- as.data.frame(signif(sig.qtable, 2))
qtable %>%
    mutate(
        `B-cell` = row.names(.),
        pvalue = cell_spec(pvalue, color = ifelse(qtable$pvalue > 0.05, "grey", "red")),
        qvalue = cell_spec(qvalue, color = ifelse(qtable$qvalue > 0.05, "grey", "red"))
    ) %>%
    kable( escape = F,
           digits = 20,
           row.names = T,
           caption = "Significant device-type q-values") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%")
Significant device-type q-values
pvalue qvalue B-cell
1 0.011 0.31 CD27-IgD+ mature naive
2 0.045 0.45 CD19+CD5+

We plotted mean levels across time for each of the B-cells that showed a statistically significant effect across devices in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the fit. Note that the reference level for the devices is the HeartMate-II (HMII).

require(reshape2, quietly = T)
df.long <- melt(df, id.vars = colnames(df)[1:13])
groups <- make.names(c("Device Type"))
names(df.long) <- make.names(names(df.long))

invisible(suppressMessages(require(Hmisc, quietly = T)))
stat_sum_df <- function(fun, geom="errorbar", ...) {
    stat_summary(fun.data = fun, geom = geom, width = 1, ...)
}

plots.ts <- mclapply(rownames(qtable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            scale_color_d3() + scale_fill_d3() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- rownames(qtable)

for(ii in c(1:length(plots.ts))){
    for(jj in 1:length(plots.ts[[ii]])){
        sumtable <- suppressMessages(summary(models.device[[rownames(qtable)[ii]]]$model))
        sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
        if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
        sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
        cat("  \n###", rownames(qtable)[ii], "\n")
        print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),,drop=F], row.names = T) %>%
                  kable_styling(bootstrap_options = c("striped", 
                                                      "hover", 
                                                      "condensed",
                                                      "responsive"),
                                font_size = 12)
        )
        cat("  \n")
        suppressWarnings(print(plots.ts[[ii]][[jj]]))
        cat("  \n")
    }
}

CD27-IgD+ mature naive

Estimate Std. Error df t value Pr(>|t|)
Device.TypePVAD -33.85486 10.06687 23.54824 -3.362998 0.0026287

CD19+CD5+

Estimate Std. Error df t value Pr(>|t|)
Device.TypePVAD 23.16366 8.615765 23.18582 2.688521 0.0130627

Linar mixed-effect model

We next analyzed the differences in B-cell levels across device types using a linear mixed effect model with time as a continuous variable, and included the interaction term. Here we report variables that had a statistically significant variance (Benjamini-Hochberg \(p<0.05\)) across devices, time, or their interaction.

suppressMessages(require(lmerTest, quietly = T))
suppressMessages(require(car, quietly = TRUE))

df.device <- df
colnames(df.device) <- make.names(colnames(df), unique = T)

varnames <- colnames(df.device)[14:42]
models.device <- mclapply(varnames, function(this_var){
    this_formula <- as.formula(paste0(this_var, " ~ Device.Type * Time + (1|PatientID)"))
    invisible(suppressMessages(this_model <- lmer(this_formula, data = df.device)))
    this_anova <- Anova(this_model, type = 2)
    pvals <- this_anova$`Pr(>Chisq)`[c(1,2,3)]
    return(list(model = this_model,
                pvals = pvals))
}, mc.cores = detectCores()-1)
names(models.device) <- colnames(df)[14:42]

pvals <- do.call(rbind, lapply(models.device, function(x) x$pvals))
rownames(pvals) <- varnames
colnames(pvals) <- c("Device", "Time", "Device:Time")

qBH <- matrix(p.adjust(pvals, method = "BH"), nrow = nrow(pvals))
rownames(qBH) <- colnames(df)[14:42]
colnames(qBH) <- colnames(pvals)

sigvars <- apply(apply(pvals, 2, function(x) x<=0.05), 1, function(x) sum(x) > 0)
sig.qtable <- qBH[sigvars,][order(apply(qBH[sigvars,], 1, min)),]

qtable <- as.data.frame(signif(sig.qtable, 2))
qtable %>%
    mutate(
        `B-cell` = row.names(.),
        Device = cell_spec(Device, color = ifelse(qtable$Device > 0.05, "grey", "red")),
        `Device:Time` = cell_spec(`Device:Time`, color = ifelse(qtable$`Device:Time` > 0.05, "grey", "red")),
        `Time` = cell_spec(`Time`, color = ifelse(qtable$`Time` > 0.05, "grey", "red"))
    ) %>%
    kable( escape = F,
           digits = 20,
           row.names = T,
           caption = "Significant device-type q-values") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%")
Significant device-type q-values
Device Time Device:Time B-cell
1 0.6 0.0023 0.007 CD19+CD11b+
2 0.5 0.025 0.82 CD268 of +27-38++transitional
3 0.82 0.025 0.82 CD19+CD5+CD11b+
4 0.85 0.089 0.82 CD19+27-38+CD5+transitionals
5 0.4 0.089 0.87 CD19+CD268+
6 0.85 0.093 0.96 CD19 of live lymph
7 0.47 0.24 0.093 CD27+38++plasma blasts
8 0.15 0.2 0.24 CD27-IgD+ mature naive
9 0.47 0.18 0.82 CD27+IgD- switched memory
10 0.26 0.19 0.93 CD19+CD5+
11 0.82 0.19 0.87 CD19+CD5+CD24hi
12 0.82 0.2 0.82 num lymph

We plotted mean levels across time for each of the B-cells that showed a statistically significant effect across devices in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit. Note that the reference level for the devices is the HeartMate-II (HMII).

require(reshape2, quietly = T)
df.long <- melt(df, id.vars = colnames(df)[1:13])
groups <- make.names(c("Device Type"))
names(df.long) <- make.names(names(df.long))

invisible(suppressMessages(require(Hmisc, quietly = T)))
stat_sum_df <- function(fun, geom="errorbar", ...) {
    stat_summary(fun.data = fun, geom = geom, width = 1, ...)
}

plots.ts <- mclapply(rownames(qtable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            scale_color_d3() + scale_fill_d3() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- rownames(qtable)

for(ii in c(1:length(plots.ts))){
    for(jj in 1:length(plots.ts[[ii]])){
        sumtable <- suppressMessages(summary(models.device[[rownames(qtable)[ii]]]$model))
        sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
        if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
        sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
        cat("  \n###", rownames(qtable)[ii], "\n")
        print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),,drop=F], row.names = T) %>%
                  kable_styling(bootstrap_options = c("striped", 
                                                      "hover", 
                                                      "condensed",
                                                      "responsive"),
                                font_size = 12)
        )
        cat("  \n")
        suppressWarnings(print(plots.ts[[ii]][[jj]]))
        cat("  \n")
    }
}

CD19+CD11b+

Estimate Std. Error df t value Pr(>|t|)
Device.TypeTAH:Time 2.1101624 0.5538816 89.91928 3.809771 0.0002541
Device.TypePVAD:Time 0.8054722 0.2943289 88.29844 2.736640 0.0075041
Time 0.2628129 0.1188074 89.80591 2.212092 0.0294973

CD268 of +27-38++transitional

Estimate Std. Error df t value Pr(>|t|)
Time -0.7873216 0.3030382 89.7994 -2.598094 0.0109556

CD19+CD5+CD11b+

Estimate Std. Error df t value Pr(>|t|)
Time 0.290629 0.1150892 91.27507 2.525249 0.0132869

CD19+27-38+CD5+transitionals

Estimate Std. Error df t value Pr(>|t|)
Time 0.0910247 0.0395736 89.94153 2.30014 0.0237538

CD19+CD268+

Estimate Std. Error df t value Pr(>|t|)
Time -0.6425061 0.230173 88.33500 -2.791405 0.0064310
Device.TypePVAD -37.7748686 16.138544 27.53028 -2.340662 0.0267323

CD19 of live lymph

Estimate Std. Error df t value Pr(>|t|)
Time -0.3683497 0.140373 93.29159 -2.624078 0.0101509

CD27+38++plasma blasts

Estimate Std. Error df t value Pr(>|t|)
Device.TypePVAD:Time 0.2768179 0.1033923 88.20169 2.677355 0.0088487

CD27-IgD+ mature naive

Estimate Std. Error df t value Pr(>|t|)
Device.TypePVAD -32.958219 10.7368535 29.72852 -3.069635 0.0045463
Device.TypeTAH:Time -2.247904 0.8166692 89.34072 -2.752527 0.0071618

CD27+IgD- switched memory

Estimate Std. Error df t value Pr(>|t|)
Time 0.3654018 0.1606941 88.62502 2.273897 0.0253875
Device.TypePVAD 25.0335542 10.8348873 28.15496 2.310458 0.0284029

CD19+CD5+

Estimate Std. Error df t value Pr(>|t|)
Device.TypePVAD 19.97858 9.433345 33.29567 2.117868 0.0417364

CD19+CD5+CD24hi

Estimate Std. Error df t value Pr(>|t|)

Two-way repeated measures anova

Finally, we attempted to analyze the differences in B-cell levels across device types using a two-way repeated measures ANOVA. Here we report variables that had a statistically significant variance (\(p<0.05\)) across devices, times, or their interaction. As there were 5 device types and 7 timepoints, but only 120 samples, this model is severely underpowered. The posterior belief in any of these results should therefore be quite small (as a consequence of Bayes rule).

suppressMessages(require(lmerTest, quietly = T))
suppressMessages(require(car, quietly = TRUE))

df.device <- df
colnames(df.device) <- make.names(colnames(df), unique = T)

varnames <- colnames(df.device)[14:42]
models.device <- mclapply(varnames, function(this_var){
    this_formula <- as.formula(paste0(this_var, " ~ Device.Type * factor(Time) + (1|PatientID)"))
    invisible(suppressMessages(this_model <- lmer(this_formula, data = df.device)))
    this_anova <- Anova(this_model, type = 2)
    pvals <- this_anova$`Pr(>Chisq)`[c(1,2,3)]
    return(list(model = this_model,
                pvals = pvals))
}, mc.cores = detectCores()-1)
names(models.device) <- colnames(df)[14:42]

pvals <- do.call(rbind, lapply(models.device, function(x) x$pvals))
rownames(pvals) <- varnames
colnames(pvals) <- c("Device", "factor(Time)", "Device:factor(Time)")

qBH <- matrix(p.adjust(pvals, method = "BH"), nrow = nrow(pvals))
rownames(qBH) <- colnames(df)[14:42]
colnames(qBH) <- colnames(pvals)

sigvars <- apply(apply(pvals, 2, function(x) x<=0.05), 1, function(x) sum(x) > 0)
sig.qtable <- qBH[sigvars,][order(apply(qBH[sigvars,], 1, min)),]

qtable <- as.data.frame(signif(sig.qtable, 2))
qtable %>%
    mutate(
        `B-cell` = row.names(.),
        Device = cell_spec(Device, color = ifelse(qtable$Device > 0.05, "grey", "red")),
        `Device:factor(Time)` = cell_spec(`Device:factor(Time)`, color = ifelse(qtable$`Device:factor(Time)` > 0.05, "grey", "red")),
        `factor(Time)` = cell_spec(`factor(Time)`, color = ifelse(qtable$`factor(Time)` > 0.05, "grey", "red"))
    ) %>%
    kable( escape = F,
           digits = 20,
           row.names = T,
           caption = "Significant device-type q-values") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 12) %>%
    scroll_box(width = "100%")
Significant device-type q-values
Device factor(Time) Device:factor(Time) B-cell
1 0.16 0.071 3e-10 CD19+CD5+
2 0.33 0.041 8.7e-06 CD27+IgD- switched memory
3 0.54 0.071 7.8e-05 CD19+CD27+
4 0.54 0.31 7.8e-05 CD27-IgD- switched memory
5 0.59 4e-04 0.0017 CD19+CD11b+
6 0.071 0.16 0.0034 CD27-IgD+ mature naive
7 0.93 0.0058 0.98 CD19 of live lymph
8 0.93 0.94 0.012 CD19CD24hiCD38-memory
9 0.27 0.041 1 CD27+IgD-IgM+ switched memory
10 0.98 0.93 0.057 CD19+24dim38dim naive mature
11 0.93 0.54 0.069 CD19+CD27+CD24hi
12 0.54 0.069 0.93 CD268 of +27-38++transitional
13 0.85 0.1 0.59 CD19+CD5+CD11b+
14 0.31 0.16 0.11 CD19+CD268+
15 0.44 0.22 0.16 CD27+38++plasma blasts

We plotted mean levels across time for each of the B-cells that showed a statistically significant effect across devices in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit. Note that the reference level for the time comparisons is timepoint 0, and the reference level for the devices is the HeartMate-II (HMII).

require(reshape2, quietly = T)
df.long <- melt(df, id.vars = colnames(df)[1:13])
groups <- make.names(c("Device Type"))
names(df.long) <- make.names(names(df.long))

invisible(suppressMessages(require(Hmisc, quietly = T)))
stat_sum_df <- function(fun, geom="errorbar", ...) {
    stat_summary(fun.data = fun, geom = geom, width = 1, ...)
}

plots.ts <- mclapply(rownames(qtable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            scale_color_d3() + scale_fill_d3() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- rownames(qtable)

for(ii in c(1:length(plots.ts))){
    for(jj in 1:length(plots.ts[[ii]])){
        sumtable <- suppressMessages(summary(models.device[[rownames(qtable)[ii]]]$model))
        sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
        if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
        sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
        cat("  \n###", rownames(qtable)[ii], "\n")
        print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),], row.names = T) %>%
                  kable_styling(bootstrap_options = c("striped", 
                                                      "hover", 
                                                      "condensed",
                                                      "responsive"),
                                font_size = 12)
        )
        cat("  \n")
        suppressWarnings(print(plots.ts[[ii]][[jj]]))
        cat("  \n")
    }
}

CD19+CD5+

Estimate Std. Error df t value Pr(>|t|)
Device.TypeTAH:factor(Time)5 79.25674 11.074689 67.92369 7.156565 0.0000000
Device.TypePVAD 35.52841 9.894557 43.53668 3.590703 0.0008326
Device.TypePVAD:factor(Time)3 -25.18076 8.751886 71.43140 -2.877181 0.0052873
Device.TypePVAD:factor(Time)5 -31.58267 11.241069 73.06105 -2.809579 0.0063611
Device.TypePVAD:factor(Time)8 -20.60508 7.572762 69.16187 -2.720946 0.0082277

CD27+IgD- switched memory

Estimate Std. Error df t value Pr(>|t|)
Device.TypePVAD:factor(Time)3 -44.11347 8.159997 70.32575 -5.406064 0.0000008
Device.TypePVAD 46.55839 10.997051 35.92281 4.233715 0.0001521
Device.TypePVAD:factor(Time)8 -25.13824 7.037695 68.87488 -3.571943 0.0006518
Device.TypePVAD:factor(Time)5 -35.76427 10.505189 71.34332 -3.404439 0.0010914
Device.TypeHVAD:factor(Time)1 22.95446 9.340788 68.84101 2.457443 0.0165160
Device.TypeTAH:factor(Time)1 19.80313 9.340788 68.84101 2.120071 0.0376068

CD19+CD27+

Estimate Std. Error df t value Pr(>|t|)
Device.TypePVAD:factor(Time)3 -44.86147 9.352169 70.81469 -4.796905 0.0000087
Device.TypePVAD 40.90474 11.677321 38.82945 3.502922 0.0011749
Device.TypePVAD:factor(Time)5 -36.38395 12.029880 72.05165 -3.024465 0.0034502
Device.TypePVAD:factor(Time)8 -22.52547 8.075489 69.06610 -2.789363 0.0068193
Device.TypeHVAD:factor(Time)1 26.28504 10.718455 69.03694 2.452316 0.0167262

CD27-IgD- switched memory

Estimate Std. Error df t value Pr(>|t|)
Device.TypeTAH:factor(Time)5 71.68679 13.81377 66.76087 5.189516 0.0000021
Device.TypeTAH:factor(Time)14 34.72518 15.34586 76.65452 2.262837 0.0264767
Device.TypePVAD:factor(Time)3 -22.85502 10.65984 75.69337 -2.144030 0.0352412

CD19+CD11b+

Estimate Std. Error df t value Pr(>|t|)
Device.TypeTAH:factor(Time)14 52.95791 10.085405 72.02616 5.250946 0.0000015
Device.TypeTAH:factor(Time)1 21.89638 8.043661 70.05677 2.722190 0.0081768
Device.TypePVAD:factor(Time)14 17.80531 7.197456 72.09536 2.473834 0.0157228
Device.TypeTAH:factor(Time)3 16.76692 8.048562 70.13474 2.083219 0.0408780

CD27-IgD+ mature naive

Estimate Std. Error df t value Pr(>|t|)
Device.TypeTAH:factor(Time)5 -49.338009 12.632153 68.31609 -3.905748 0.0002179
Device.TypePVAD -40.588774 11.742845 41.95252 -3.456469 0.0012664
Device.TypeTAH:factor(Time)14 -42.518754 14.393799 71.18856 -2.953963 0.0042493
Device.TypeTAH:factor(Time)1 -29.588950 11.470396 69.40423 -2.579593 0.0120122
factor(Time)8 -7.421838 2.996063 68.79529 -2.477197 0.0157031
Device.TypeHVAD:factor(Time)1 -23.895494 11.470396 69.40423 -2.083232 0.0409155

CD19 of live lymph

Estimate Std. Error df t value Pr(>|t|)
factor(Time)1 8.670607 2.838543 69.69728 3.054598 0.0031913
Device.TypeHVAD:factor(Time)5 18.621453 8.106642 68.64256 2.297061 0.0246707
Device.TypeHVAD:factor(Time)3 16.313101 7.998551 68.60746 2.039507 0.0452529

CD19CD24hiCD38-memory

Estimate Std. Error df t value Pr(>|t|)
Device.TypeTAH:factor(Time)5 -22.909686 6.066396 68.15378 -3.776490 0.0003363
Device.TypeHVAD:factor(Time)1 12.942595 5.519995 68.69011 2.344675 0.0219407
factor(Time)1 -3.388861 1.568014 68.34505 -2.161244 0.0341830

CD27+IgD-IgM+ switched memory

Estimate Std. Error df t value Pr(>|t|)
factor(Time)5 10.711338 3.308708 71.05832 3.237317 0.0018346
factor(Time)3 8.467757 3.003200 71.42351 2.819578 0.0062189
factor(Time)14 8.273763 3.705650 70.67012 2.232742 0.0287354

CD19+24dim38dim naive mature

Estimate Std. Error df t value Pr(>|t|)
Device.TypePVAD:factor(Time)3 -51.79949 13.95422 76.12859 -3.712101 0.0003892

CD19+CD27+CD24hi

Estimate Std. Error df t value Pr(>|t|)
Device.TypeTAH:factor(Time)5 -22.54139 9.215036 68.32836 -2.446154 0.0170182
factor(Time)8 4.88642 2.184243 68.95972 2.237123 0.0285135
Device.TypeHVAD:factor(Time)1 18.31135 8.355636 69.79036 2.191497 0.0317555

CD268 of +27-38++transitional

Estimate Std. Error df t value Pr(>|t|)
factor(Time)8 -13.21005 5.828370 69.05529 -2.266508 0.0265563
factor(Time)14 -16.08481 7.874549 69.65369 -2.042632 0.0448754

CD19+CD5+CD11b+

Estimate Std. Error df t value Pr(>|t|)

CD19+CD268+

Estimate Std. Error df t value Pr(>|t|)
Device.TypeTAH:factor(Time)5 -66.88010 17.16044 67.82057 -3.897342 0.0002253
Device.TypePVAD -39.35919 17.41639 37.46865 -2.259893 0.0297274

CD27+38++plasma blasts

Estimate Std. Error df t value Pr(>|t|)
Device.TypeTAH:factor(Time)1 7.049114 2.979632 70.86462 2.365767 0.0207319
Device.TypeCMAG 7.332777 3.630995 44.46774 2.019495 0.0494876

HeartMate-II analysis

The HeartMate-II (HMII) recipients were the largest group, and we analyzed them by themselves due to the previously observed variability across devices.

require(reshape2, quietly = T)
df.HMII <- subset(df, df$`Device Type`=="HMII")
df.long <- melt(df.HMII, id.vars = colnames(df)[1:13])

groups <- make.names(c("AgeGreater60", 
                       "Sex",
                       "LowIntermacs",
                       "RVAD", 
                       "Sensitized",
                       "VAD Indication", 
                       #"Device Type", 
                       "Survival",
                       "Outcome"))

names(df.long) <- make.names(names(df.long))

plots.ts <- mclapply(unique(df.long$variable), function(this_var){
    lapply(groups, function(this_groups){
        ggplot(subset(df.long, df.long$variable == this_var)) +
            aes(x = Time, y = value, group = PatientID) +
            aes_string(color = this_groups, fill = this_groups) +
            geom_line(alpha = 0) + 
            geom_point(alpha = 0) + 
            #stat_smooth(aes_string(group = this_groups), method = "loess", span = 1) +
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("point"), position = position_dodge(.5)) + 
            stat_summary(fun.y = mean, aes_string(group = this_groups), geom=c("line"), position = position_dodge(.5)) + 
            stat_sum_df(function(x) mean_cl_normal(x, conf.int = 0.68), mapping = aes_string(group = this_groups), position = position_dodge(.5)) + 
            scale_color_aaas() + scale_fill_aaas() +
            xlab("Time (days after surgery)") +
            ylab(this_var) + 
            ggtitle(paste(this_var)) +
            theme_classic()
    })
}, mc.cores = detectCores()-1)
names(plots.ts) <- unique(df.long$variable)

# for(ii in c(1:length(plots.ts))){
#     for(jj in 1:length(plots.ts[[ii]])){
#         suppressWarnings(print(plots.ts[[ii]][[jj]]))
#     }
# } 

PCA

Using PCA, we found large variability between individual patients, compared to the variability within individual patients. None of the other features were clearly separable.

suppressMessages(require(ggbiplot, quietly = T))
require(ggsci, quietly = T)

# Efron's double standardization
double_standardize <- function(x, niter = 100) {
    for(i in 1:niter) x <- t(scale(t(scale(x))))
    return(as.data.frame(x))
}

isna <- unique(unlist(apply(df.HMII[,14:42], 2, function(x) which(is.na(x)))))
pca <- prcomp(double_standardize(df.HMII[-isna, 14:42]), center = TRUE, scale. = TRUE)

colorfun <- function(grouping, ...){
    if(nlevels(factor(grouping)) > 10) scale_color_discrete(...)
    if(nlevels(factor(grouping)) <=10) scale_color_d3(...)
} 

plots.pca <- mclapply(names(df.HMII)[1:17], function(this_var){
    this_groups <- df.HMII[-isna, this_var]
    if(is.numeric(this_groups)) this_groups <- cut(this_groups, 4)
    
    ggbiplot(pca,
             groups = this_groups,
             ellipse = TRUE,
             alpha = 0.3,
             varname.size = 1.2) + 
        colorfun(this_groups, name = this_var) +
        ggtitle(paste0(this_var, " variability")) +
        theme_classic() 
}, mc.cores = detectCores()-1)

names(plots.pca) <- names(df.HMII)[1:17]
#plots.pca$PatientID
#plots.pca$`Device Type`

for(ii in 1:length(plots.pca)){
    cat("  \n###", names(plots.pca)[ii], "\n")
    suppressWarnings(print(plots.pca[[ii]]))
    cat("  \n")
}

PatientID

Age

AgeGreater60

Sex

LowIntermacs

InterMACS

RVAD

Sensitized

VAD Indication

Device Type

Survival

Outcome

Time

num Total PBMC

num lymph

lymph

live lymph

One-way repeated measures anova

We analyzed the differences in B-cell levels for various features using a one-way repeated measures anova. Here we report variables that had a statistically significant variance (\(p<0.05\)) across groups, or groups at each timepoint.

suppressMessages(require(lmerTest, quietly = TRUE))
suppressMessages(require(car, quietly = TRUE))
require(reshape2, quietly = TRUE)

df.lmer <- df.HMII
names(df.lmer) <- make.names(names(df.lmer), unique = TRUE)

groupvars.ix <- c(3,4,5,7,8,9,11)
groupvars <- names(df.lmer)[groupvars.ix]

bcells.ix <- c(14:42)
bcells <- names(df.lmer)[bcells.ix]

models.b <- mclapply(groupvars, function(this_groupvar){
    models.bcells <- lapply(bcells, function(this_bcell){
        this_formula <- as.formula(paste0(this_bcell, " ~ ", this_groupvar, 
                                          " + (1|PatientID)"))
        suppressMessages(suppressWarnings(this_model <- lmer(this_formula, data = droplevels(df.lmer))))
        this_anova <- Anova(this_model, type = 2)
        this_pvalues <- this_anova$`Pr(>Chisq)`
        names(this_pvalues) <- rownames(this_anova)
        #return(this_pvalues)
        return(list(model = this_model,
                    pvals = this_pvalues))
    })
    names(models.bcells) <- colnames(df)[14:42]
    pvalues <- do.call(rbind, lapply(models.bcells, function(x) x$pvals))
    rownames(pvalues) <- bcells
    #return(pvalues)
    return(list(model = models.bcells,
                pvals = pvalues))
}, mc.cores = detectCores()-1)
names(models.b) <- groupvars
pvals <- lapply(models.b, function(x) x$pvals)
# something wrong here
names(pvals) <- groupvars
pvals.matrix <- do.call(cbind, lapply(pvals, function(this_pval) this_pval[,c(1)]))



# Benjamini Hochberg
# qBH <- matrix(p.adjust(as.numeric(pvals.matrix), 
#                        method = "BH"), 
#               nrow = nrow(pvals.matrix), 
#               ncol = ncol(pvals.matrix), 
#               byrow = F) 
# rownames(qBH) <- rownames(pvals.matrix)
# colnames(qBH) <- colnames(pvals.matrix)
# rownames(qBH) <- names(df)[bcells.ix]
# qvalsBH.df <- melt(qBH)
# colnames(qvalsBH.df) <- c("B-cell", "parameter", "qvalue")
# qvalsBH.df.ranked <- qvalsBH.df[order(qvalsBH.df$qvalue, decreasing = F),]
# qvalsBH.df.ranked[qvalsBH.df.ranked$qvalue <= 0.3,]

# Local FDR
require(fdrtool, quietly = T)
invisible(suppressMessages(fdrobj <- fdrtool(as.numeric(pvals.matrix), statistic = "pvalue", plot = F, verbose = F)))
qvals.matrix <- matrix(fdrobj$q, nrow = nrow(pvals.matrix), ncol = ncol(pvals.matrix), byrow = F)
rownames(qvals.matrix) <- rownames(pvals.matrix)
colnames(qvals.matrix) <- colnames(pvals.matrix)
rownames(qvals.matrix) <- names(df)[bcells.ix]

pvals.df <- melt(pvals.matrix)
qvals.df <- melt(qvals.matrix)
colnames(qvals.df) <- colnames(pvals.df) <- c("B-cell", "parameter", "qvalue")
qvals.df.short <- qvals.df[pvals.df$qvalue <= 0.05,]
shortlist <- qvals.df.short[order(qvals.df.short$qvalue),]



kable(shortlist, 
      digits = 3,
      row.names = T,
      caption = "Significant results") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 10) %>%
    scroll_box(width = "100%")
Significant results
B-cell parameter qvalue
87 CD19+27+IgD-38++IgG ASC LowIntermacs 0.104
32 lymph Sex 0.314
2 num lymph AgeGreater60 0.348
3 lymph AgeGreater60 0.417
45 CD27+IgD+IgM+ nonswitched memory Sex 0.675
101 CD27+IgD+ unswitched memory RVAD 0.721
29 CD19+27+IgD-38++IgG ASC AgeGreater60 0.725
20 CD19CD24hiCD38-memory AgeGreater60 0.732
61 lymph LowIntermacs 0.732
55 CD19+CD27+CD24hi Sex 0.738

We plotted the average across time for each of the B-cells that showed a statistically significant effect across various factors in the above one-way anova. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit.

require(stringr, quietly = T)
siggroups <- sapply(str_split(shortlist$parameter, ":"), function(x) x[1])
for(ii in 1:nrow(shortlist)){
    this_group <-siggroups[ii]
    this_bcell <- as.character(shortlist$`B-cell`[ii])
    cat("  \n###", as.character(shortlist$`B-cell`[ii]), "\n")
    
    sumtable <- suppressMessages(summary(models.b[[this_group]]$model[[this_bcell]]$model))
    sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
    if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
    sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
    
    print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),,drop=F], row.names = T) %>%
              kable_styling(bootstrap_options = c("striped", 
                                                  "hover", 
                                                  "condensed",
                                                  "responsive"),
                            font_size = 12)
    )
    cat("  \n")
    
    suppressWarnings(print(plots.ts[[shortlist$`B-cell`[ii]]][[which(groups == this_group)]]))
    cat("  \n")
}

CD19+27+IgD-38++IgG ASC

Estimate Std. Error df t value Pr(>|t|)
LowIntermacsHigh 2.078718 0.60011 17.90374 3.463894 0.0027878

lymph

Estimate Std. Error df t value Pr(>|t|)
SexMale -18.44662 6.370603 17.92503 -2.895584 0.0096664

num lymph

lymph

Estimate Std. Error df t value Pr(>|t|)
AgeGreater60older -15.22801 5.790664 17.53477 -2.629751 0.0172536

CD27+IgD+IgM+ nonswitched memory

Estimate Std. Error df t value Pr(>|t|)
SexMale -16.26667 7.449772 17.76457 -2.183512 0.0426655

CD27+IgD+ unswitched memory

Estimate Std. Error df t value Pr(>|t|)

CD19+27+IgD-38++IgG ASC

Estimate Std. Error df t value Pr(>|t|)

CD19CD24hiCD38-memory

Estimate Std. Error df t value Pr(>|t|)

lymph

Estimate Std. Error df t value Pr(>|t|)

CD19+CD27+CD24hi

Estimate Std. Error df t value Pr(>|t|)

Linear mixed-effect model

We analyzed the differences in B-cell levels for various features using a linear mixed effect model, with time as a continuous variable. Here we report variables that had a statistically significant variance (\(p<0.05\)) across groups, or groups at each timepoint.

suppressMessages(require(lmerTest, quietly = TRUE))
suppressMessages(require(car, quietly = TRUE))
require(reshape2, quietly = TRUE)

df.lmer <- df.HMII
names(df.lmer) <- make.names(names(df.lmer), unique = TRUE)

groupvars.ix <- c(3,4,5,7,8,9,11)
groupvars <- names(df.lmer)[groupvars.ix]

bcells.ix <- c(14:42)
bcells <- names(df.lmer)[bcells.ix]

models.b <- mclapply(groupvars, function(this_groupvar){
    models.bcells <- lapply(bcells, function(this_bcell){
        this_formula <- as.formula(paste0(this_bcell, " ~ ", this_groupvar, 
                                          " * Time + (1|PatientID)"))
        suppressMessages(suppressWarnings(this_model <- lmer(this_formula, data = droplevels(df.lmer))))
        this_anova <- Anova(this_model, type = 2)
        this_pvalues <- this_anova$`Pr(>Chisq)`
        names(this_pvalues) <- rownames(this_anova)
        #return(this_pvalues)
        return(list(model = this_model,
                    pvals = this_pvalues))
    })
    names(models.bcells) <- colnames(df)[14:42]
    pvalues <- do.call(rbind, lapply(models.bcells, function(x) x$pvals))
    rownames(pvalues) <- bcells
    #return(pvalues)
    return(list(model = models.bcells,
                pvals = pvalues))
}, mc.cores = detectCores()-1)
names(models.b) <- groupvars
pvals <- lapply(models.b, function(x) x$pvals)
# something wrong here
names(pvals) <- groupvars
pvals.matrix <- do.call(cbind, lapply(pvals, function(this_pval) this_pval[,c(1,3)]))



# Benjamini Hochberg
# qBH <- matrix(p.adjust(as.numeric(pvals.matrix), 
#                        method = "BH"), 
#               nrow = nrow(pvals.matrix), 
#               ncol = ncol(pvals.matrix), 
#               byrow = F) 
# rownames(qBH) <- rownames(pvals.matrix)
# colnames(qBH) <- colnames(pvals.matrix)
# rownames(qBH) <- names(df)[bcells.ix]
# qvalsBH.df <- melt(qBH)
# colnames(qvalsBH.df) <- c("B-cell", "parameter", "qvalue")
# qvalsBH.df.ranked <- qvalsBH.df[order(qvalsBH.df$qvalue, decreasing = F),]
# qvalsBH.df.ranked[qvalsBH.df.ranked$qvalue <= 0.3,]

# Local FDR
require(fdrtool, quietly = T)
invisible(suppressMessages(fdrobj <- fdrtool(as.numeric(pvals.matrix), statistic = "pvalue", plot = F, verbose = F)))
qvals.matrix <- matrix(fdrobj$q, nrow = nrow(pvals.matrix), ncol = ncol(pvals.matrix), byrow = F)
rownames(qvals.matrix) <- rownames(pvals.matrix)
colnames(qvals.matrix) <- colnames(pvals.matrix)
rownames(qvals.matrix) <- names(df)[bcells.ix]

pvals.df <- melt(pvals.matrix)
qvals.df <- melt(qvals.matrix)
colnames(qvals.df) <- colnames(pvals.df) <- c("B-cell", "parameter", "qvalue")
qvals.df.short <- qvals.df[pvals.df$qvalue <= 0.05,]
shortlist <- qvals.df.short[order(qvals.df.short$qvalue),]



kable(shortlist, 
      digits = 3,
      row.names = T,
      caption = "Significant results") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 10) %>%
    scroll_box(width = "100%")
Significant results
B-cell parameter qvalue
31 num lymph AgeGreater60:Time 0.041
145 CD19+27+IgD-38++IgG ASC LowIntermacs 0.041
2 num lymph AgeGreater60 0.056
347 CD19+CD5+CD11b+ VAD.Indication:Time 0.060
105 CD19+24dim38dim naive mature Sex:Time 0.217
61 lymph Sex 0.247
344 CD19+CD5+ VAD.Indication:Time 0.319
155 CD27-38++ transitional LowIntermacs:Time 0.331
343 CD19+CD11b+ VAD.Indication:Time 0.358
3 lymph AgeGreater60 0.401
346 CD19+CD5+CD24hi VAD.Indication:Time 0.413
53 CD19+CD11b+ AgeGreater60:Time 0.439
30 num Total PBMC AgeGreater60:Time 0.508
39 CD27-38++ transitional AgeGreater60:Time 0.537
382 CD3 of live lymph Survival:Time 0.549
74 CD27+IgD+IgM+ nonswitched memory Sex 0.549
92 CD3 of live lymph Sex:Time 0.569
97 CD27-38++ transitional Sex:Time 0.572
188 CD27+IgD+ unswitched memory RVAD 0.574
119 lymph LowIntermacs 0.594
57 CD19+CD5+CD11b+ AgeGreater60:Time 0.605
20 CD19CD24hiCD38-memory AgeGreater60 0.613
29 CD19+27+IgD-38++IgG ASC AgeGreater60 0.615
328 CD27+38++plasma blasts VAD.Indication:Time 0.616
84 CD19+CD27+CD24hi Sex 0.618
321 num lymph VAD.Indication:Time 0.622
52 CD268 of +27-38++transitional AgeGreater60:Time 0.626
399 CD19+CD268+ Survival:Time 0.627

We plotted the average across time for each of the B-cells that showed a statistically significant effect across various factors in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit.

require(stringr, quietly = T)
siggroups <- sapply(str_split(shortlist$parameter, ":"), function(x) x[1])
for(ii in 1:nrow(shortlist)){
    this_group <-siggroups[ii]
    this_bcell <- as.character(shortlist$`B-cell`[ii])
    cat("  \n###", as.character(shortlist$`B-cell`[ii]), "\n")
    
    sumtable <- suppressMessages(summary(models.b[[this_group]]$model[[this_bcell]]$model))
    sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
    if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
    sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
    
    print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),,drop=F], row.names = T) %>%
              kable_styling(bootstrap_options = c("striped", 
                                                  "hover", 
                                                  "condensed",
                                                  "responsive"),
                            font_size = 12)
    )
    cat("  \n")
    
    suppressWarnings(print(plots.ts[[shortlist$`B-cell`[ii]]][[which(groups == this_group)]]))
    cat("  \n")
}

num lymph

Estimate Std. Error df t value Pr(>|t|)
Time 8206.115 1791.316 1287.688 4.581054 0.0000051
AgeGreater60older:Time -8800.327 2381.169 2274.895 -3.695802 0.0002244

CD19+27+IgD-38++IgG ASC

Estimate Std. Error df t value Pr(>|t|)
LowIntermacsHigh 3.058696 0.774233 45.35361 3.950614 0.0002697

num lymph

Estimate Std. Error df t value Pr(>|t|)
Time 8206.115 1791.316 1287.688 4.581054 0.0000051
AgeGreater60older:Time -8800.327 2381.169 2274.895 -3.695802 0.0002244

CD19+CD5+CD11b+

Estimate Std. Error df t value Pr(>|t|)
VAD.IndicationDT:Time 0.7535101 0.2208725 66.53183 3.411516 0.0011029

CD19+24dim38dim naive mature

Estimate Std. Error df t value Pr(>|t|)
Time -0.8161430 0.2628953 64.83737 -3.104441 0.0028253
SexMale:Time 0.9045098 0.3076620 64.86978 2.939946 0.0045445

lymph

Estimate Std. Error df t value Pr(>|t|)
SexMale -17.61531 7.156886 25.50523 -2.46131 0.0209393

CD19+CD5+

Estimate Std. Error df t value Pr(>|t|)
VAD.IndicationDT:Time 0.7495301 0.2761314 67.17233 2.714397 0.0084299

CD27-38++ transitional

Estimate Std. Error df t value Pr(>|t|)
Time 0.2728189 0.0943080 68.33055 2.892849 0.0051174
LowIntermacsHigh:Time -0.3029055 0.1126997 67.21515 -2.687721 0.0090598

CD19+CD11b+

Estimate Std. Error df t value Pr(>|t|)
VAD.IndicationDT:Time 0.5911116 0.2251629 65.39624 2.625262 0.0107693

lymph

Estimate Std. Error df t value Pr(>|t|)
Time 0.7428399 0.3610038 69.6636 2.057706 0.0433624

CD19+CD5+CD24hi

Estimate Std. Error df t value Pr(>|t|)
VAD.IndicationDT:Time 0.4514841 0.1801168 68.30036 2.506618 0.0145744

CD19+CD11b+

Estimate Std. Error df t value Pr(>|t|)
Time 0.5499941 0.1571744 65.94444 3.499260 0.0008412
AgeGreater60older:Time -0.5051776 0.2059164 65.50982 -2.453314 0.0168239

num Total PBMC

Estimate Std. Error df t value Pr(>|t|)
Time 8539.807 3421.824 8916.676 2.495688 0.0125892
AgeGreater60older:Time -10420.156 4502.010 15842.034 -2.314556 0.0206499

CD27-38++ transitional

Estimate Std. Error df t value Pr(>|t|)
Time 0.1924463 0.0792502 69.27151 2.428340 0.0177687
AgeGreater60older:Time -0.2343863 0.1043092 67.73025 -2.247034 0.0279003

CD3 of live lymph

Estimate Std. Error df t value Pr(>|t|)
Survivaldead:Time -1.829251 0.8242817 66.32952 -2.219206 0.0298979

CD27+IgD+IgM+ nonswitched memory

Estimate Std. Error df t value Pr(>|t|)
SexMale -20.7260608 7.6448089 21.84155 -2.711129 0.0128054
Time -0.9006538 0.3895893 66.29131 -2.311803 0.0239045

CD3 of live lymph

Estimate Std. Error df t value Pr(>|t|)
SexMale:Time -1.143545 0.5277063 66.36504 -2.167011 0.0338279

CD27-38++ transitional

Estimate Std. Error df t value Pr(>|t|)
Time 0.2418057 0.0994538 65.99621 2.431336 0.0177663
SexMale:Time -0.2512294 0.1163638 66.05080 -2.159000 0.0344873

CD27+IgD+ unswitched memory

Estimate Std. Error df t value Pr(>|t|)
RVADYes 8.138537 3.426992 21.41969 2.374834 0.0269627

lymph

Estimate Std. Error df t value Pr(>|t|)

CD19+CD5+CD11b+

Estimate Std. Error df t value Pr(>|t|)
Time 0.5382593 0.1602283 67.49846 3.359327 0.0012883
AgeGreater60older:Time -0.4334330 0.2101801 66.78306 -2.062198 0.0430827

CD19CD24hiCD38-memory

Estimate Std. Error df t value Pr(>|t|)

CD19+27+IgD-38++IgG ASC

Estimate Std. Error df t value Pr(>|t|)
AgeGreater60older 2.044292 0.8970807 39.13157 2.278828 0.0282137

CD27+38++plasma blasts

Estimate Std. Error df t value Pr(>|t|)
Time 0.1018402 0.0463711 65.98620 2.196200 0.0315970
VAD.IndicationDT:Time -0.1781890 0.0880343 66.32797 -2.024086 0.0469921

CD19+CD27+CD24hi

Estimate Std. Error df t value Pr(>|t|)

num lymph

CD268 of +27-38++transitional

Estimate Std. Error df t value Pr(>|t|)

CD19+CD268+

Estimate Std. Error df t value Pr(>|t|)
Time -0.5331745 0.2039893 64.53338 -2.613738 0.0111329

Two-way repeated measures anova

We analyzed the differences in B-cell levels for various features using a two-way repeated measures ANOVA. Here we report variables that had a statistically significant variance (\(p<0.05\)) across groups, or groups at each timepoint.

suppressMessages(require(lmerTest, quietly = TRUE))
suppressMessages(require(car, quietly = TRUE))
require(reshape2, quietly = TRUE)

df.lmer <- df.HMII
names(df.lmer) <- make.names(names(df.lmer), unique = TRUE)

groupvars.ix <- c(3,4,5,7,8,9,11)
groupvars <- names(df.lmer)[groupvars.ix]

bcells.ix <- c(14:42)
bcells <- names(df.lmer)[bcells.ix]

models.b <- mclapply(groupvars, function(this_groupvar){
    models.bcells <- lapply(bcells, function(this_bcell){
        this_formula <- as.formula(paste0(this_bcell, " ~ ", this_groupvar, 
                                          " * factor(Time) + (1|PatientID)"))
        suppressMessages(suppressWarnings(this_model <- lmer(this_formula, data = droplevels(df.lmer))))
        this_anova <- Anova(this_model, type = 2)
        this_pvalues <- this_anova$`Pr(>Chisq)`
        names(this_pvalues) <- rownames(this_anova)
        #return(this_pvalues)
        return(list(model = this_model,
                    pvals = this_pvalues))
    })
    names(models.bcells) <- colnames(df)[14:42]
    pvalues <- do.call(rbind, lapply(models.bcells, function(x) x$pvals))
    rownames(pvalues) <- bcells
    #return(pvalues)
    return(list(model = models.bcells,
                pvals = pvalues))
}, mc.cores = detectCores()-1)
names(models.b) <- groupvars
pvals <- lapply(models.b, function(x) x$pvals)
# something wrong here
names(pvals) <- groupvars
pvals.matrix <- do.call(cbind, lapply(pvals, function(this_pval) this_pval[,c(1,3)]))



# Benjamini Hochberg
# qBH <- matrix(p.adjust(as.numeric(pvals.matrix), 
#                        method = "BH"), 
#               nrow = nrow(pvals.matrix), 
#               ncol = ncol(pvals.matrix), 
#               byrow = F) 
# rownames(qBH) <- rownames(pvals.matrix)
# colnames(qBH) <- colnames(pvals.matrix)
# rownames(qBH) <- names(df)[bcells.ix]
# qvalsBH.df <- melt(qBH)
# colnames(qvalsBH.df) <- c("B-cell", "parameter", "qvalue")
# qvalsBH.df.ranked <- qvalsBH.df[order(qvalsBH.df$qvalue, decreasing = F),]
# qvalsBH.df.ranked[qvalsBH.df.ranked$qvalue <= 0.3,]

# Local FDR
require(fdrtool, quietly = T)
invisible(suppressMessages(fdrobj <- fdrtool(as.numeric(pvals.matrix), statistic = "pvalue", plot = F, verbose = F)))
qvals.matrix <- matrix(fdrobj$q, nrow = nrow(pvals.matrix), ncol = ncol(pvals.matrix), byrow = F)
rownames(qvals.matrix) <- rownames(pvals.matrix)
colnames(qvals.matrix) <- colnames(pvals.matrix)
rownames(qvals.matrix) <- names(df)[bcells.ix]

pvals.df <- melt(pvals.matrix)
qvals.df <- melt(qvals.matrix)
colnames(qvals.df) <- colnames(pvals.df) <- c("B-cell", "parameter", "qvalue")
qvals.df.short <- qvals.df[pvals.df$qvalue <= 0.05,]
shortlist <- qvals.df.short[order(qvals.df.short$qvalue),]



kable(shortlist, 
      digits = 3,
      row.names = T,
      caption = "Significant results") %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive"),
                  font_size = 10) %>%
    scroll_box(width = "100%")
Significant results
B-cell parameter qvalue
2 num lymph AgeGreater60 0.040
145 CD19+27+IgD-38++IgG ASC LowIntermacs 0.040
52 CD268 of +27-38++transitional AgeGreater60:factor(Time) 0.073
276 CD27+IgD-IgM+ switched memory Sensitized:factor(Time) 0.081
61 lymph Sex 0.200
348 CD19+27+IgD-38++IgG ASC VAD.Indication:factor(Time) 0.213
232 CD19+27+IgD-38++IgG ASC RVAD:factor(Time) 0.230
155 CD27-38++ transitional LowIntermacs:factor(Time) 0.254
3 lymph AgeGreater60 0.257
347 CD19+CD5+CD11b+ VAD.Indication:factor(Time) 0.377
31 num lymph AgeGreater60:factor(Time) 0.378
101 CD27+IgD+ unswitched memory Sex:factor(Time) 0.417
40 CD27-IgD+ mature naive AgeGreater60:factor(Time) 0.438
119 lymph LowIntermacs 0.528
321 num lymph VAD.Indication:factor(Time) 0.592
74 CD27+IgD+IgM+ nonswitched memory Sex 0.616
217 CD27+IgD+ unswitched memory RVAD:factor(Time) 0.621
103 CD27+IgD+IgM+ nonswitched memory Sex:factor(Time) 0.632
225 CD19+CD268+ RVAD:factor(Time) 0.643
399 CD19+CD268+ Survival:factor(Time) 0.675
344 CD19+CD5+ VAD.Indication:factor(Time) 0.677
224 CD19+27-38+CD5+transitionals RVAD:factor(Time) 0.705
51 CD19+CD268+ AgeGreater60:factor(Time) 0.705
188 CD27+IgD+ unswitched memory RVAD 0.708
29 CD19+27+IgD-38++IgG ASC AgeGreater60 0.727
20 CD19CD24hiCD38-memory AgeGreater60 0.732

We plotted the average across time for each of the B-cells that showed a statistically significant effect across various factors in the above mixed effect models. We drew attention to specific features that induced the positive test result, by listing the model parameters with \(p<0.05\) in the multivariate fit.

require(stringr, quietly = T)
siggroups <- sapply(str_split(shortlist$parameter, ":"), function(x) x[1])
for(ii in 1:nrow(shortlist)){
    this_group <-siggroups[ii]
    this_bcell <- as.character(shortlist$`B-cell`[ii])
    cat("  \n###", as.character(shortlist$`B-cell`[ii]), "\n")
    
    sumtable <- suppressMessages(summary(models.b[[this_group]]$model[[this_bcell]]$model))
    sumtable <- as.data.frame(sumtable$coefficients)[-1, ,drop=F] # drop intercept
    if(!("Pr(>|t|)" %in% colnames(sumtable))) next()
    sigsum <- sumtable[sumtable[,"Pr(>|t|)"] <= 0.05, , drop = F]
    
    print(kable(sigsum[order(sigsum[,"Pr(>|t|)"]),,drop=F], row.names = T) %>%
              kable_styling(bootstrap_options = c("striped", 
                                                  "hover", 
                                                  "condensed",
                                                  "responsive"),
                            font_size = 12)
    )
    cat("  \n")
    
    suppressWarnings(print(plots.ts[[shortlist$`B-cell`[ii]]][[which(groups == this_group)]]))
    cat("  \n")
}

num lymph

Estimate Std. Error df t value Pr(>|t|)
factor(Time)21 169798.9 44963.12 1780.755 3.776403 0.0001643
AgeGreater60older:factor(Time)21 -214015.6 62340.16 1651.364 -3.433029 0.0006116
AgeGreater60older:factor(Time)5 -136885.7 61168.57 2679.738 -2.237844 0.0253131
factor(Time)14 128214.1 57636.41 1470.433 2.224532 0.0262639

CD19+27+IgD-38++IgG ASC

Estimate Std. Error df t value Pr(>|t|)
LowIntermacsHigh 3.095 1.279721 69 2.418497 0.0182277

CD268 of +27-38++transitional

Estimate Std. Error df t value Pr(>|t|)
AgeGreater60older:factor(Time)8 -27.22177 10.33554 54.58048 -2.633803 0.0109633
factor(Time)14 -24.92941 12.43679 55.44507 -2.004489 0.0499111

CD27+IgD-IgM+ switched memory

Estimate Std. Error df t value Pr(>|t|)
SensitizedYes:factor(Time)3 17.965877 5.302290 23.31974 3.388324 0.0024949
factor(Time)5 9.747852 4.547781 23.30369 2.143431 0.0427391

lymph

Estimate Std. Error df t value Pr(>|t|)
SexMale -22.22807 9.001586 53.16972 -2.46935 0.0167811

CD19+27+IgD-38++IgG ASC

Estimate Std. Error df t value Pr(>|t|)
factor(Time)5 4.722429 1.330913 60.14685 3.548263 0.0007592
VAD.IndicationDT:factor(Time)5 -6.577445 1.874482 57.09611 -3.508940 0.0008850

CD19+27+IgD-38++IgG ASC

Estimate Std. Error df t value Pr(>|t|)
RVADYes:factor(Time)5 11.55470 2.912778 59.01744 3.966900 0.0001997
factor(Time)14 -2.17979 1.043488 56.60858 -2.088945 0.0412202

CD27-38++ transitional

Estimate Std. Error df t value Pr(>|t|)
LowIntermacsHigh:factor(Time)21 -11.380113 2.894440 59.65377 -3.931715 0.0002222
factor(Time)21 9.440964 2.450218 61.33982 3.853112 0.0002815

lymph

Estimate Std. Error df t value Pr(>|t|)

CD19+CD5+CD11b+

Estimate Std. Error df t value Pr(>|t|)
VAD.IndicationDT:factor(Time)21 18.60843 5.699112 56.68547 3.265145 0.0018584
VAD.IndicationDT:factor(Time)8 10.32734 4.156816 55.37116 2.484436 0.0160275

num lymph

Estimate Std. Error df t value Pr(>|t|)
factor(Time)21 169798.9 44963.12 1780.755 3.776403 0.0001643
AgeGreater60older:factor(Time)21 -214015.6 62340.16 1651.364 -3.433029 0.0006116
AgeGreater60older:factor(Time)5 -136885.7 61168.57 2679.738 -2.237844 0.0253131
factor(Time)14 128214.1 57636.41 1470.433 2.224532 0.0262639

CD27+IgD+ unswitched memory

Estimate Std. Error df t value Pr(>|t|)
SexMale:factor(Time)8 9.848232 2.746704 54.41204 3.585472 0.0007200
factor(Time)8 -8.159681 2.352215 54.37377 -3.468936 0.0010300
factor(Time)14 -8.101354 3.062814 55.12449 -2.645069 0.0106199
SexMale -10.353852 3.922084 26.14988 -2.639885 0.0138005
factor(Time)5 -7.781354 3.062814 55.12449 -2.540590 0.0139123
SexMale:factor(Time)5 8.682914 3.494502 55.06965 2.484736 0.0160327
SexMale:factor(Time)3 5.898797 2.849382 54.68689 2.070202 0.0431688
SexMale:factor(Time)14 7.448953 3.625289 54.98976 2.054720 0.0446698

CD27-IgD+ mature naive

Estimate Std. Error df t value Pr(>|t|)
factor(Time)14 -19.34727 6.403688 54.77954 -3.021270 0.0038210
AgeGreater60older:factor(Time)14 17.50253 7.722534 54.66849 2.266424 0.0274049

lymph

Estimate Std. Error df t value Pr(>|t|)

num lymph

CD27+IgD+IgM+ nonswitched memory

Estimate Std. Error df t value Pr(>|t|)
SexMale -31.50752 9.007301 37.48741 -3.497997 0.0012251
SexMale:factor(Time)8 27.02760 8.735159 54.73057 3.094116 0.0031071
SexMale:factor(Time)3 24.44017 9.046112 55.34852 2.701732 0.0091386
factor(Time)8 -20.00903 7.482252 54.65698 -2.674199 0.0098595
factor(Time)3 -15.59653 7.482252 54.65698 -2.084470 0.0418040

CD27+IgD+ unswitched memory

Estimate Std. Error df t value Pr(>|t|)
RVADYes:factor(Time)8 -10.16212 3.121996 55.38516 -3.255008 0.0019363
RVADYes 11.57901 3.999725 34.08235 2.894952 0.0065730

CD27+IgD+IgM+ nonswitched memory

Estimate Std. Error df t value Pr(>|t|)
SexMale -31.50752 9.007301 37.48741 -3.497997 0.0012251
SexMale:factor(Time)8 27.02760 8.735159 54.73057 3.094116 0.0031071
SexMale:factor(Time)3 24.44017 9.046112 55.34852 2.701732 0.0091386
factor(Time)8 -20.00903 7.482252 54.65698 -2.674199 0.0098595
factor(Time)3 -15.59653 7.482252 54.65698 -2.084470 0.0418040

CD19+CD268+

Estimate Std. Error df t value Pr(>|t|)
RVADYes:factor(Time)3 -27.99543 10.480918 56.66602 -2.671086 0.0098514
RVADYes:factor(Time)8 -20.33650 9.017817 55.39973 -2.255146 0.0280950

CD19+CD268+

Estimate Std. Error df t value Pr(>|t|)
Survivaldead:factor(Time)8 -21.73092 8.069861 55.34955 -2.692849 0.0093560
Survivaldead:factor(Time)3 -18.76413 8.701780 55.67137 -2.156355 0.0353912

CD19+CD5+

Estimate Std. Error df t value Pr(>|t|)
VAD.IndicationDT:factor(Time)8 13.77699 5.156496 55.63351 2.671773 0.0098785
VAD.IndicationDT:factor(Time)21 18.57915 7.056032 57.34507 2.633087 0.0108559

CD19+27-38+CD5+transitionals

Estimate Std. Error df t value Pr(>|t|)

CD19+CD268+

Estimate Std. Error df t value Pr(>|t|)
factor(Time)14 -34.03871 8.813153 54.50989 -3.862263 0.0003004
AgeGreater60older:factor(Time)14 34.33335 10.627446 54.41429 3.230631 0.0020977
factor(Time)21 -17.19694 6.912037 54.73282 -2.487970 0.0159221
AgeGreater60older:factor(Time)3 19.32183 8.006711 54.60648 2.413204 0.0191971
factor(Time)8 -11.89000 5.511085 53.89561 -2.157470 0.0354450

CD27+IgD+ unswitched memory

Estimate Std. Error df t value Pr(>|t|)
RVADYes:factor(Time)8 -10.16212 3.121996 55.38516 -3.255008 0.0019363
RVADYes 11.57901 3.999725 34.08235 2.894952 0.0065730

CD19+27+IgD-38++IgG ASC

Estimate Std. Error df t value Pr(>|t|)

CD19CD24hiCD38-memory

Estimate Std. Error df t value Pr(>|t|)